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АВ5ТКАСТ 


A method for the computation of confidence intervals for 
Circular error probability (CEP) based on first order variance 
estimates was introduced in 1966. It was later found that 
under certain conditions the resulting confidence intervals 
for CEP were smaller than expected. As a result a second 
order variance estimate method was developed, at the Johns 
Hopkins University Applied Physics Laboratory, which greatly 
improved the accuracy of the confidence intervals for CEP. 
The purpose of this thesis is to develop and test procedures 
for the 3-dimensional case to obtain a second order estimate 


for variance of spherical error probability (SEP). 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed 
in this research may not have been exercised for all cases of 
interest. While every effort has been made, within the time 
available, to ensure that the programs are free of 
computational and logic errors, they cannot be considered 
validated. Any application of these programs without 


additional verification is at the risk of the user. 
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I. INTRODUCTION 


A. BACKGROUND 

In 1966 W.R. Blischke and A.H. Halpin (Ref. 1] first 
described methods for the computation of confidence intervals 
for circular error probability (CEP) based on first order 
variance estimates. When these procedures were implemented at 
the Johns Hopkins University Applied Physics Laboratory (JHU- 
APL) it was found that under certain conditions the resulting 
confidence intervals for CEP were smaller than expected. As 
a result R.C. Ferguson, K.V. Kitzman and P.B. Jackson [Ref. 2] 
have extended the Blischke-Halpin procedures to create a 
second order variance estimate. Testing conducted by Kitzman 
[Ref. 3] has confirmed that this method produces a more 
accurate estimate for the variance of CEP. 

Our goal will be to extend this procedure to the 3- 
dimensional case to obtain a second order estimate for 


variance of spherical error probability (SEP). 


В. METHODOLOGY 

We begin by assuming that missile detonations are 
distributed as trivariate Gaussian. In the past, analysis of 
ballistic missile test firing data has failed to disprove this 


assumption. If the mean and variance are given by: 


Hi 


Нз 


01: 61: Vaa 
2 = [01 022 9255]. 


03: G32 Oaa 


Then the probability of a detonation occurring within a 


distance R of the origin is given by: 


АК; р, >) = I и 
D 


where, 


X1 
x = X, , 
X 
а? = аец], 


and D - (x,,2,,2x,)|xi*x2*x3 s R?). 


By the implicit function theorem, < > O implies that the 


equation AR;p,2) = CONSTANT defines a differentiable function 


Кр, Ў) such that AR(p,2);p, 2) = CONSTANT. If the CONSTANT is 


set to 2. then 2) 1s the SEP function. Let X,,X,,..,X, be 


a random sample of independent and identically distributed 


random vectors with distribution N(u,2). Then the maximum 
1 n 
likelihood estimators of ы апа X аге шары: апа 
i= 
1 п 
- = eG :-0). When fi and X replace y and X in the 
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equation Z(R;yp,X) - then R(fi,2) gives SEP, an estimator of 


SEP. 
The second-order approximation for the variance of SÊP is 


given by: 


OSep = (DyyuSEP)7P(D,yuSEP) 


+ <р ЗЕР) (PB (DA. SEP) 


-Holsen s] a 


1709 Т ° ° ° 
where X" 7 (0,,,0,,,0,,,022,,02,,0,,) , D,, y.SEP is the derivative of 


SEP with respect to н," viewed as a column vector, D'SuSEP is 


the usual second derivative matrix (or Hessian) of SEP with 


PUES 
respect ton n E - eee is the variance-covariance 
Pipe Pou 


matrix R, 5", & is the tensor product operation and the 


underline notation is used to represent the column vector 


formed by stringing out the elements of the matrix by rows. 
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For example, if A= 
pie, 3 4 


then A -[12 3 4]T. Note that fi is 








distributed N(u, ŠE) END dn =. It can also be seen thatn2 
has a Wishart distribution with parameter X, so that 
Py = = (282) . The complete derivation of equation (1), based 


on Taylor Series expansion, is given in appendix A. 


Since ñ and Š are independent [Ref. 4:p. 102] it follows 


ПАС Р, «= 0 and equation (1) can be written as: 


E (ЕР) e (seen. eye { 2522) 
SEP = _ди_ p Qu. apu > XE 














21 | see) E E | eve ў р. ® | E 
cr — sue JU EE one 











92чи. -Oz"Qu) duon" -Quon" 





| ELS ES (232) ер, NE 











&sEP|? 925ЕР)/ (2) 
AEE o) e) 


Derivation of the first order terms in equation (2) has 
been accomplished by S.D. Hill at JHU-APL [Ref. 5]. In order 


to solve for the second order terms it will be necessary to 


find expressions for F SEP E о ЕЕ апа SEP i.e. it is 
ди? дӯ ч Ən o>" 











necessary to батар, СЕР, the second derivative matrix of SEP 


with respect to its arguments. Chapter II will describe the 
derivation of the second order variance estimate and Chapter 
III will discuss implementation and testing of the algorithm. 


Conclusions will be presented in Chapter IV. 


ІІ. DERIVATION OF THE SECOND ORDER ESTIMATE 


A. FIRST ORDER TERMS 
We begin by describing the derivation of the first order terms 


given in equation (2). If we define: 


Ңх,,х;,хұ;р, X") = ехрі-> (ж-н) 15 arn) 


where 


сії 012 gi? 
L-1 = |621 022 ç23| , (13 


631 632 633 


l H I 
а 
W 
к 


(X-u) °Ч(х-н)= (х,-н,)%011+2(х,-н,)(х„-н›)0°%+2(х, -нь,)(х,-н;)о3* 
+(x, оо +2(х, G -Из)672 +(x3 TL 033 
апа 


д(р,0,ф;ы, 27) = Ервіпбсовф,рвіпбвеіпф,рсов0Ө;ы,239. 


Then, with a shift to spherical coordinates, we see that: 


AR; p, D4) = — 1 — f. [, f, e*sin9(p,0,6: v, Z9apad8a8. 


(/2n)'q (3) 


We wish to solve equation (3) for each of the partials 


OSEP 
where А = (Hy, Mor B31 F447 F427 9437 922: 023-033) - From the 


OA ; 


1 





R, 2) = is the SEP function, so we need 


previous section, Б 


to solve for the partials = 





i 
Taking partials of equation (3) we see that, 


дю __ 1 OF 
ae ae = Creel lag (4) 


and 


PR) - TRE "віпӨо(к,Ө,ф)аӨсф. 
Gaara |. 


We can solve the right hand side of equation (4) for 


i=1,2,3 by using the relationships 


99 =- 99 (99... 99 99 --.94 
ди, Ox, Ор; Ox, ди, ох, 








зо ме see that 


-J[[ 32 &aamax, = [[[ 32 pa aao 


If we apply Green's Theorem to this equation and then make 


a shift to spherical coordinates we get 


Л; og = dx, dx, = R?["*[”cosbsin'0gR, 0, ф)аВаф 


f [[ 58 ax ax ax, = R*[™ [*sinosin?0gqR, 8, >)d8dd 
D 


SN Ax, dx, dx, = 5? |" ["созбзз лв, 0, $) ава 


z °" [*sin20ə(R, Ө, ф)аӨаф 


То solve the right hand side of equation (4) for i=4 + 39 


we first note that: 








2-26 ‚ dej 
C E 
Op,0p, dg N 1%7 
00 55 


зо ме get 








у ОВОС 
о , i= 
2 dx; Op; y 
д да pn. 
I š Jl ao 7) 
OX; 9р; 


which is used to evaluate the integral of og as follows: 


ff [E a aas 


Jil ox dx,dx,dx, 
D 

- ГГГ-99 ах. ах. ах, 
ПГ 


a 09. ax, dx, dx, 
ШЕ; 


-J[f 38 tas 





On, 


1 


2 л гл 
= N f, [o (psin0cosó-p,)*o"*(psin0sin$-y;) 


*o^(pcos0-p,)|cosósin*6g(n, Ө,ф)а9аф 


2x nmn N . ° 
gi f [o**(psinOcoso-p,)+o**(psinOsind-p,) 


+0°*(pcos6-p,)|sindsin*6 GR, 6, >) dbdb 


2 N fn 
Ef J: (6-Қ(рвіпӨсозвф-н,)ғ0 “((рвіпӨвзіпф-н.) 


+o**(pcosO-p,)|sin26q(R, 0, ф) 9 аф 


2 к гл 
RU f [o?*(psin6cos-p,)+o7*(psinOsind-p,) 
0 0 


%0%(рсов0-н,)|віпфвзіп”Өо(в,0,%ф)адаф 


2 N (K 
A [Tor psindcont-ayjePpsindeiag-ny 
0 


*o^*(pcos0-p,)|sin20s(n, 6, фуабаф 


Е д9 а R* 2* nn N W 32 : . = 
50 аа, 2 f. JA |67(рвіпӨсовф H,)+o°*(psinOsing-p,) 


*o**(pcos0-p,)|sin20g(R, 6, ф)а8аф 


If we define the functions: 


CC, MUS. D ГР Геов!00%сов51Ф)4Е,8,ф)404ф 
CS, j,k (Ri, D") = |” |Үсов()6)зіп1%)4Е,0,ф)484% 
SCi, j,k (Rip, D”) = [ "| віп()0)совҚ1)4(Е,8,ф)а84ф 


SS, k (Ri 2) = [ "| sin'G6sin*u4)g, 0, ф) аф 


Extending the work of Blischke and Halpin, with the aid of 


Green's Theorem, we can see that: 


OSERp, DY) _ | 56,1,1,1 552,1,1,1 SC1,2,0,0 à 


(5) 
Op 5С1,1,0,0 ЭС1,1,0,0 25С1,1,0,0 (ey, 29); y, 29) 
[1 г В, 
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95С, 1,1,1 
др, 

955, 11,1 
Ән, 

OSC, , 0,0 


OSERp, DY) _ 1 a 


(6) 
one 2 SC, 1,0,0 055, 111 


ди; 
д5С, 20,0 

ды. 
OSC, 2,0,0 

20p, 


2 


(eu ruya, zt). 


Applying the partials given in equation (6) to the 


partials -[[[ SZ ax, ax,ax, given above we see that: 
D 1 


OSEP ШЕГЕ 11 12 13/6, = 
Soe онанизм, w) 


- и. (о*1 р, +012 ,+013,)} 





90, К 


OSEP _ А) 


мос р, +0124, +013) 


OSEP — 1 


до,» 2 W; 





(2 07 (5 79) *o (s 9) eos ws 


-w,(o?p, *ot?p, +013 ,}} 
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ОБЕРЕ 22 23. 
n. ume M 


-Wio(0!? p, *0??u,*0*?p,)) 











OSEP т | к 

Ser "oma (zie m ma) + om 0) +070, +] 
wo, +0? p, +073, )) 

OSEP ere 

Ser "az [o c) mona 


-w,(o23p, +073, +073p,)} 


where WR;p,2") is given by: 


W СС, 1,1,1 

м - СС, 3,1,1 

wa = С5 1,1,1 

W, = CS1,3,1,1 

Ws = SC1,1,0,0 

We = SC1,2,0,0 
R Ва 

7 1,370 O 

Ws = SC;,1i,1,1 

Ws 7 SC, 12,1 

Ито = 55, 1,1,1 

Мі 7 55, 1,1,2 


Мі2 ә 55, 12,1 


where each of the functions is evaluated at Мы,2%,р, 25. 


This completes the derivation of the first order terms 


required in the computation of the variance estimate. 
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В. CONCEPT FOR DERIVATION OF THE SECOND ORDER TERMS 


We first observe that equations (5) and (6) enable one to 


express ps as a composition of three functions 


во. ап, each Of which is a fairly simple function of its 


arguments. 


Explicitly, 


D, swSERp, ZY) = T,°T,°T (hp, E") 


where, 


H 
ТҰРУЫ Е 2 
SERV, 2") 
H 
yu 
Т.һ, 24, R) = > 
WR; p, 2”) 
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Wg 
















W5 
mo 
Ws 
We 
2 Ws 
W 
H Е (ow 5204 (иаа) а ыш АЫ 
5 
2 R dollw .+4012w..+03(w.-w б cilu., +gł2u. +01? 
n'|-| x09 isse (ката тала ШІ LN 
5 
W R M 
qu (9 (s -w)*o t (ws но ил - 2 (ор, +о?р.+о??р;) 
5 5 
R dä 
S (2075 W11 +405 W t0 (w,7w,))- 29 и Г. 
5 5 
K 


W 
(c (и, А +о22(0, -w,)*o?*(w, +Ws)) E ll р, +022 H2 +923 H3) 





4 W; 
к міз, оз о MEE LM ME 
av (9 (WTW) +O (wo и) +07 (w; ws) (o to^ uz tou) 
5 5 | 


Since D, ySEAp, 2") = T,°T,°T,(u, 2") we may apply the chain 


rule to obtain the second derivative as a matrix product: 


Da ,SSEHMu,ZI") = DI(Z TIKE) KDE ERI ee Е (7) 


The simplicity of the functions Т,, 7, and 7, means that 
the derivatives DT,, DT, and DT, are easily obtained and then 
the rather large matrix products involved in DT,DT;DT, can be 
left to the computer leading to a rather painless computation 


2 
of Ы, ұ«СЕР. 
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ШІ (еһе аевіуаєіме of the function, T,) is the 10x9 


matrix consisting of the 9x9 identity matrix (from Б 
p, i" 


with the tenth row consisting of the partials of SEP with 


respect to mu and sigma. 


DT, 


О ешо O OROM 
© © >). [ei 7. © 
О ОО CLE O 
о ооо P оО 
о эро 
ооо ооо 
ооо оо 
Fx O O O O O O O 
о оо оооло о 


0 0 0 0 0 0 0 0 
ОСЕР дЅЕР дЅЕР ОСЕР ОСЕР ОСЕР ОСЕР ОСЕР ОСЕР 


Op, дн; ды, 00, 90,, да; до, до, 005; | 


= 


DT, is the 22x10 matrix consisting of the 10x10 identity 


и 
matrix (from Í . , with rows 11 to 22 consisting of the 
H, , 


12х10 sub-matrix of partials of W with respect to ы, 24 апа 


R. 
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Dit, = 


о о о о о о о о о H 

e o o o oro Cfo kc 

5. 
о о о о о о о нр o o 
о о о о о ор о о о 
о о о о о рг о о о о 
о о о о гю о о о о o 
о о о гю о о о о о o 
о орг о о о о о о о 
ою о о о о о о о o 
е о о о о о о о о o 








DT, is the 9x22 matrix consisting of the partials of 


u;(i71,..,9) with respect to y, 2", R and W. 


а 1х3 ыш. 1х6 I 1x1 E Ix12 
T ) 3g ) Sp ) SW X12) 


DT. 
ди, 


ox" 





2 1x3 1x6 SE 1x1 JE 1x12 
ОШ ) (1х6) Sr | ) "ay ) 


Thus, in order to obtain the second degree variance 
estimate given in equation (2) we must solve for the values of 


all elements of DI D and DIE Once these values are 


computed, the solution is found by the matrix multiplication 
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shown in equation (7), yielding the 9x9 matrix of second 
partials of SEP with respect to mu and sigma. 


The values of the DT, elements are computed from equations 








(5) and (6). The determination of DT, requires the evaluation 
of 58, D. and = which is fully described in appendix B. 


Finally, to solve for the elements of DT, it is necessary to 


solve for each of the partials given above. The evaluation of 


these partials is given in appendix C. 
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III. IMPLEMENTATION AND TESTING 


А. PROGRAMMING 


The calculation of o%», the second order approximation for 


variance of SEP (when р апа Ў аге known) has been implemented 
in the FORTRAN program SEPCMP which is provided in appendix D. 
The following steps are followed to obtain the approximation: 


- Input the distribution mean (H) and variance- 


covariance matrix elements (2%. 

- Compute SEP using H and 2. 

- Compute the values of each element of the matrices 
DT1, DT2 and DT3 and then form the matrix product 
comprising the second derivative of SEP with 


respect to p and XZ". The first derivative of SEP 


with respect to р апа Ў" іѕ then given by row 10 
of the matrix DT1. 
- Calculate the variance-covariance matrices of fi 


(P) and of E (Py) from 2. 
- Form the SEP variance approximation using equation 


(2). 

Computation of SEP given p and 2"! is based on the PL-1 
program RAP provided by JHU-APL. RAP is a general program 
which provides the radius of coverage for any probability in 
the interval (0,1) for either 2 or 3 dimensions. Those 
portions of the program applicable to this problem, ie. 
probability - 0.5 and 3 dimensions, have been used in the 
FORTRAN subroutine FINDSEP to obtain SEP.  Derivation of the 
equation used to compute SEP is given by L. S. Simpkins [Ref. 


6] and is implemented in the subroutine FINDSEP using Gaussian 
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Quadrature (via the IMSL subroutine QTWODQ) to evaluate the 
required double integral. 

Gaussian Quadrature is also used to evaluate each of the 
double integrals in the W, CC, CS, SC and SS functions, which 
are required in the evaluation of the matrices DT1, DT2 and 


DT3. 


B. TESTING 

The purpose of deriving the second order approximation of 
the variance of SEP is to yield more accurate confidence 
intervals for SEP. The approximation of these confidence 


intervals is derived from the fact that since fi and X are 


maximum likelihood estimators of p and X and thus SÉP is 


asymptotically Normal with mean SEP. Thus, for example, an 


approximate 95% confidence interval for SEP is given by 


SËP+1 ..96 /var(SËP). 


In order to test the affect of using the second order 
approximation for the variance of SEP, tests of confidence 
interval coverage were conducted as follows: 


- Sample 30 random vectors from Np,2) and compute 
û апа Ê for the sample. 


- Calculate SEP from ñ and 2. 


- Calculate 6%,, by replacing fi and Í for р апах 


in both the first order and second order 


approximations for the variance of SEP. 
- Compute the large sample approximation 95% 


confidence intervals for SEP using SÉP£1.96/625. 
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- Repeat 1000 times to calculate the percentage of 
times that confidence intervals based on each 
approximation cover the true SEP. 

- For each input variance-covariance matrix, D, 
repeat the test at different magnitudes of the 
mean. 

Trivariate normal random vectors are obtained using the 
IMSL subroutines DCHFAC and DRNMVN. DCHFAC performs Cholesky 
factorization of the distribution SIGMA which is then used as 
an input to DRNMVN which returns the desired random vector. 


The estimates of P, and Py are found by replacing H and 
У with f апа Ê in the formulations of P, and P, given in 


Chapter I. Тһе results of this test are summarized in Table 
ШІ 


TABLE 1. 95 PCT CONFIDENCE INTERVAL COVERAGE 


| SIGMA MU FIRST ORDER | SECOND ORDER 
| CI CI 
5.2 


















10,000 3,000 2,000 
3,000 10,000 3,000 
2,000 3,000 10,000 
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These results indicate that there is little accuracy to be 
gained by using the second order variance approximation. 
There is, however, a situation in which the second order 
approximation becomes important. This case has been described 
by K. V. Kitzman [Ref. 3] and is paraphrased here. 

Weapons system targeting can be measured from separate 
subsystems (such as navigation, guidance, postboost, etc.). 
Suppose that there are n independent subsystems and let 

ла th 


Х:1,....Хш = Җи, Ў) represent the errors for the i subsystem 


i=1,...,n, then the impact errors for the weapons system are 


n e e 
y Ja 

ze = Жі; 1-1, еее ,m and W ec. J у = Мру=пр, 201122) o 
1=1 


If the Y's are measured directly and used to estimatep, 


20 n 
and 2, (to get the system SEP) then в. = = 1 = and 


2 
Pa = £D Ej = “1 (EF). If on the other hand the individual 
wn у т 


n 
X's are available from the subsystems then Й, = Yy X; where 
ta 


п 

v 1 ү A A п 

K. = = ж X;; and 2. š = 3 (X:5 BY X ü). Then B, - =» and 
321 


Zu 


a 22082). The significant factor is that the estimate of 


Py, differs from the first estimate by a factor of n while the 
estimate of Puy remains unchanged. Thus as the number of 


subsystems becomes larger and the number of observations 
becomes smaller the size of P, can become very small іп 


relation tor As the mean approaches O it follows that 


ОЗЕР 
ды 





approaches 0 as well and thus the contribution of P, to 


the estimate can become insignificant even though it is large 


in relation to the P, term. 


This affect has been examined as follows: 


- Divide р апа ХУ ру п, the number of iid subsystems 
being simulated. 

- Draw m random samples from the distribution given 
by this new u and £ for each of the n subsystems. 

- Calculate fi and Ê for each of the n subsystems. 


A 


- Sum the n estimates to get a total fi and Ê, the 
estimate for mean and variance in Y. 


- Compute confidence intervals for SEP uSing both 
the first and second order approximations with 


ТА _ 2n 
Pi = Ro and Py, = — (282). 


This test was conducted with n = 4 апа т = 5. Тһе results 
of this test are given in Table 2. As can be seen, the effect 
is greatest when the mean is 0, and decreases as the mean 


becomes larger, which is what would be expected. 
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TABLE 2. 95 PCT CI COVERAGE, MULTIPLE SUBSYSTEMS 


SIGMA MU FIRST ORDER | SECOND ORDER 
ЕТ | 
| 100 -160 270 92.4 
-160 400 -480 
270 -480 900 
5 91.4 


10 
20 
















To test the sensitivity of the second order variance 


approximation procedure to the assumption that the detonations 
are normally distributed, tests were conducted in which random 
vectors were drawn from a mixed normal distribution. The test 
is conducted in the following manner: 
- Prior to each sample obtain random variable p from 
U(0,1) distribution. 
-— If p is less than .98 then draw a sample from 
N(u,2) and proceed. 
- if p is greater than .98 then draw a sample from 
Np,10*2) and proceed. 
- Continue until 30 samples have been drawn, then 
compute all values and proceed as above. 
This procedure simulates a distribution with the same mean 
but larger tails than the normal distribution. The confidence 
interval coverages are then compared to those obtained in the 


base case from the same input distribution. The results of 


this test are summarized in Table 3. 
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TABLE 3. 95 PCT CI COVERAGE, MIXED DISTRIBUTION 


FIRST ORDER | SECOND ORDER 
CI 


10,000 3,000 2,000 
3,000 10,000 3,000 
2,000 3,000 10,000 
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IV. CONCLUSIONS 


Examination of the results in Table 1 shows that the 
second order approximation procedure provides a slight 
improvement in the 95% confidence interval coverage. The real 
utility of this method is brought out by the results of the 
multiple subsystem test shown in Table 2, where confidence 
interval coverage improves by as much as 5%. Since this case 
represents the actual situation in missile testing (i.e. very 
small mean and independent system observations) it is 
reasonable to state that use of this procedure will improve 
the estimate of the variance in SEP. Furthermore it seems 
highly unlikely that there would be any advantage to be gained 
by pursuing a higher order estimate of the variance of SEP. 

The mixed distribution test shows that the second order 
approximation is only slightly more robust than the first 
order approximation if the distribution is actually from a 
mixed and not a true normal. It also indicates that if the 
variance is extremely large then the affect of the mixing is 
greater than for the smaller values and that only a few 
extreme values are required to greatly reduce the confidence 


interval for the estimate. 
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APPENDIX A 


In order to derive the equation for the second order 
approximation for variance of SEP (equation (1)) we start by 
noting that the second order Taylor Series Expression for the 


SEP function is given by (ignoring the remainder term): 


SERA) = SERA)+D, p SERÀ) TAA + = AAD? pSERA)AA 


where A, D, ySEP, and D; ySEP are defined in the text, 4 is the 


vector consisting of the maximum 1ikelihood estimators for 


each element of A and ДА = А-А. Our goal then is to find an 


approximation for the variance of SEP(À). Because SEP(Ä) is 


asymptotically unbiased and the mean square error 


E(SER(A)-SER))?]| » Var(SER(A)) «(E SER(A)]- SER). , 


we can get this approximation by looking at E[SEP(À)-SEP(A)]?. 
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Applying the tensor identity ABC = AWC7B to the last term 


and then transposing we can write the second order Taylor 


equation as: 


SEKR)-SERA) = [D, :5ЕҢАУЛТАА» (A ATBAAN[D? z SEP) 


= [Dy ,SEP(A)]7AA+= [Dz 8ЕҢ2) ҚАХФА2) 


Now we note that since AA) =A it follows that ДАДА) = 0. 


If we square each side and take the expected value we see 


that: 


ESERA)-SERA)| * [D, y«SER)] "HA XA AT[D, yuSERA)] 


d 


+= [DX г ЗЕҢ) ТЕ(ААФАХХАХФА №", zS ER) 


In order to evaluate H(AAWAA)AAWAA)") we note that 
var(AA@Ad) = 2(SR) (POP) (SR)™ [Ref. 7:p. 44] where the matrix 


SR has the following qualities (Ref. 7: pp. 32-38]: 
- SR is symmetric 
- (SR) (SR) = SR 
- 2(SR)(u@y = u®v+ vu 
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for example, if u and v are of dimension 2 then 


140NNO. 0 
ЕО 
282 
SR = A 
EL ^0 
2 2 
о ni 


Then since, 


E(AXC9A AA AC9A A)7] 2 var(A ACA A)* ELA AC9 A AYE(A AC9A A] 7 


and 


ДАДАТ = AJOA 


if we define P to be HAAAAT)- var(AÀA) we see that 


var(AACOAA)*E(AACOA AE(A AC9A A)]T = 2 SR(POP)SR)T+Pp PT 


2(SR)'(EGS9P(SR)*P PT. 
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Since Df sSEP is a symmetric matrix it follows that 


(SR)D, pSEP - D; ;SEP and we see that: 


HSERA)-SERA)| * [D, y SER(A)]T HD, s SER()] 
-2Іш Е5ЕҢА)| POHDI 5ЕҢ2) 
+ [р SERA)" PPÍD? ЗЕҢ], 


which in turn is approximately equal to the variance of 
SEP(À). 

Finally, since X is a 3X3 symmetric matrix, we need only 
describe the unique terms (which we have denoted by X") to 


arrive at the second order equation for SEP given by equation 


(1). 
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АРРЕМОТХ В 


А. APPROACH 


It is seen that each of the 12 W, functions consists of 


a double integral of cos and sin functions multiplied by the 
function g(Rsinbcosd,Rsinösind,Rcos0). Each of the required 
derivatives carries through the integrals, so ee we need to 
evaluate each in respect to the function g. 

We will first evaluate each of the partials of g, and then 


carry these through to the functions W,. 


Recall, 


(R, 0,6) = exp|- 2|(Rsin6cos-u)*a** 


+ 


2(Rsin8cos-p,)\(Rsin6singd-p,)o** 


ai 


2(Rsin0coso-p,y(Rcos0-p,)o*? 
* (Rsin0sinó$-g;yo** 

+ 2(RSin0sin$-p;(Rcos0-p,)o*? 
* (Kcos0-y,)*o??]). 
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If we define the vectors A,, A,, A, and A, by: 


sin?0cos?$ 
2віп2бсовфвіпф 
— 2sinOcos8cosd 
_ sin20sin2p |’ 
2sin0cos0sino| 
cos?b 


ѕіпдсоѕфи, 
sin8singdp, +sinOcosdoyp, 
cos0p, *sinO0cosQu, 
Аг = 510051пфн, 
cos6p,+sinOsindgp, 


соѕӨџ, 
A, = [Wi 2р,р, 28,6, и 2n;u, uj], 


RsinOcoso-p, 
E Rsin0sin$-u,!.. 
RcosO-p, 
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Then treating R as an independent parameter and taking the 


partial with respect to R we have: 


og TEE i -usi 
= = -94{0:{РзіпӨсоз?ф-р,зіпдсозф) 


%01%(2Евіп”ӨӘсовфвіпф-нвіпӨвіпф-н,віпӨсовф) 


%015(2ЕвіпӨсовӨсовф-н,сов0-н;віпӨсовф) 


| 

( 
+о22(реіп^Өѕіпф-р,ѕіпдѕіпф) 
%0%% 2 Евіпбвіпфсов0-н;сов0-н;51п0віпф) 
| 


%033(Есов%0-н;со50)) 


= (2) [R:4,-4)]. 


Similarly, we derive: 





og . g(o**(Rsin0coso-p,) “01% ЕвіпӨвіпф -H2) +0'*(Rcosb -H3)} 


= gla 612 013] д} 


roc = g{o1*(Rsin8coso-p,)+077(RsinOsind-p,)+07>(Rcos6-p,)} 


Op, 
= 901? 9? o: j) 
25. - g (o *(Rsin8cos$-p,)*o?(Rsin0sin$-p;)*o?*(Rcos0-p,)) 
| 


т go? g23 o33]A,) 
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In order to derive the partials with respect to sigma we 


досі] 
90,; 


will first generate the partials 





and the partials of = 


with respect to 0,,. These values are given below. 





















































gil 012 013 022 
G 
д go nu - 01291! NOS Ес: 
95,1 a 
-0; G 
=o oo = 933 7512912 23 794613512 E2012 
G 
12 а q 
G -0 220 
5° TOUS = "lS Еца n 13 SB Ш> 222613 
G 
13 q q q 
-0 
56 a Peg 2: Ша. = ECT, 22022 
G 
22 q q 
=2а G G 
002, a g? qg- 
G G G 
x Е E ligi x Boos E1334 u- В: 
O33 а а а 
023 033 2-1 
= 1 
д 923 _ 23,11 922 33,11 Ou 
"=p С 039 
90, q q 24 
G -20 -6012 
E 2 23512 zu - E5312 
9015 q q q 
E 
- 2923913 -2033013 6 
0014 q q 
G -g?? 
д -g?3g?? Tn -033022 
00, а? 24 
E Жоға 
11 -2g?3g?? -2g?3g?? 
00; а? q 
д -g?? 
e 23933 Bos 
003; 2q 
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We can now use these values to find the desired partials: 


og 


дс 





11 


*(Reos8-p,)}{ 00-222] 


= 9: 2(бвіпбсовф-ц,)(о:3)4 


+2(RsinOcos-p,)\(RsinOsing-p Ja g 


) 
*2(Rsin8cosó-p,(Rcos8-pn Оа 


*(Rsin0sin$-p;) it 022011- өз) 


О 
*2(Rsin0sin$-p;(Rcos0 егемен) 
q 





a 


git 
2 


g{ 5 [04 R(2-4) "A, -2Ro* (D4) “A, +o) “A, ] 


sfo o o 


d: 


E 
q 


cmo oo ctm 9% 
2 


G G G 


iE а (к 2А,-2КА, A] 


2 G g? 


G 


= 11 
ел, ока, az 


б а 


2 


q q 


Following this same methodology it is easy to see that: 








G -0 -0 G 
longano — — 9 — зі pena, 22s ao 


Dos 204 20 а 








2 


2q° 29. а = 


-0 G с _ 
oo EO =) (о 23 22 13 = o] [x 2A, -2 RA, Ae) 


ak 
2 


i 


-0 G -0 22 
loa 33 0 13 00 бы рау санлы = | 


а“ Ge 
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29 gjon u, 923 біз 912 g бор [R2A, -2RA, +A,]-07° 
go... ЕЖ 2а2 2а? 20° 


-0 G -g 3з 
— 22 2 o —u oo [R2A, -2RA,+A,] — 
q^ а? q* 2 


dg ү 
0933 г 





g33(X 1) “| 





We are now able to calculate each of the partials of thew, 


by simply computing the effects of the four A vectors on the 
original W functions and applying the formulas derived on the 
previous pages. We will carry through the complete derivation 


for W, and then show the solutions for the other cases. Since 
the vector A, does not involve trigonometric functions, it 


will not be necessary to compute any partials for it. 


B5 


В. PARTIALS OF W1 


W = СС, т ГГ Геовбсовфжв,8,0)28аф во ме сап see that 


the integrals of gA will yield: 


sin?0cosO0cos?$ 
2sin?0cosOcos?$sino 
rp UN Е ША 215% + 
o Jo о Јо sin?0cosO0cos$sin?$ 
| 2sinOcos?0sin$cos$ 
сов3бсовф 
СС, 1,3,1 СС, 1,3,1 
2(С5, 11,17С65;,1,1,1 7 С51,1,3,15СӚ; 1,3,1) 
2(5С1,1,2,1—5С3,1,2,1) 
СС, 1,1,1 СО;,1,1,1 7 СС,1,3,17 СС; 1,3,1 
55 1,1,2 95 1,1,2 
СС, 1,1,1 


зіпӨсоѕӨсоз?фуи, 
ѕзіпдсоѕӣсоѕфѕзіпфуы, +5іпдсоѕдсоѕ?фуџ, 
соѕ?Өсоѕфџ, +5іпдсоѕдсоз? 
ГГ созбсовфоға, - 05 | Hs | Фи», 
о 70 o Jo sinbcosbcosdsindyu, 
cos?Bcosbu,+sinbcosbcosdsindu, 


соѕ?Өсоѕфи, 


2.56, 2,2,1 
H155,,2,1,2* 2M25C1,2,2,1 
48. СС, 1,1,1+2135С;,2,2,1 

№255; 2,1,2 
455CC, 1,1,1*M355,,2,1,2 

4р3СС, 1,1,1 


= 2 
4 
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ЮвіпбсовӨбсов%ф-совбсовфи, 
ГГ совбсозфғға, = Ws Rsin0cos0cos$sin$-cos0cosóy, 
0 0 0 0 


Ксов“Өсовф-совӨсовфи; 


15іп20сов?ф-р,созсозф 


—sin20sin20-n,cos0cos$ 
Рсов%бсовф-н;совбсовф 
Ін. 


a 
E $51,2,1,2R- 485 CC 4, 1,1, 
4аСС, 11,415-4р3СС1,1,1 


For future use, we point out that the partials of A, andA, 


involve only 3 distinct trigonometric terms, which we will 


define by: 


1 = sin6cosd 
B = (8, = ѕіпӨѕіпф 
В, = соѕӨ 


then we see that 


Bip, 
Boh, +p, 
Вр. +В. В 

В.р, 
Bp; * Bua 

В, 
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апа 


Thus, we need only compute the effects of A, and B on each 


of the 12 W functions to obtain the desired partials. 


C. PARTIALS OF W2 


2n ru 
W, — CC. NE J. J. сов3Өсовфо(Е,0,ф)абаф so the partials are 


given by: 


віп2бсов3бсов3ф 
2sin?0cos30cos?d$sinó 
2sin0cos0cos30cos?o 
sin?0cos30sin?$coso 
2віпбсовбсов3бвіпфсовф 
cos?0cos30coso 


cos30cosQ$A, - 


2СС, ҙ,з,1 СС 5,3,1 СС 1,3,1 
2(2С5, ; 11-С515,1,1:С51,11,172С5, ҙ,;,1%С5, ,3,1%С51 13.1) 
2(5С, 5,2,17561,1,2,1) 
2СС 3,1,1 CCa,5,1,1 CC1,1,1,12CC1,3,3,1 + CC1,5,3,1* CC1,1,3.1 
5251,5,1,2 991,1,1,2 
СС, 5,1,172 СС, 3,1,1 + CC1,1,1.1 


E 
4 
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віпбсов3бсов2ф 2(5С, 4,2,1 7 5С,,2,2,1) 
сов3бсовфВ - |віпбсов3бсовфвіпфі - i 55, 4,1,279512,1,2 
совбсов3бсовф 2(СЕ 88 2, .) 


D. PARTIALS OF W3 


2х рх 
Wee eS) TA И j^ cos0sin$sg(R,0,$)d0db so the partials are 


given by: 


sin?0cosO0cos?$sin$ 

2віп2бсовбсовфвіп?ф 

. 2віпбсов2бсовфвіпф 

совӨвіпфА, - | f 

sin20cos0sin°$ 

2sinO0cos?0sin?$ 
cos?0sinó$ 


CS11,1,1 CS3,1,1,1 CS1,1,3,1 *CS3,1,3,1 
2(CC; 1,1,37CC3,1,3,17 CC5,1,1,1* C€3,1,5,1) 
55,,1,1,27953,1,1,2 
С5, 1,3,17С5ҙ,1,3,1 
2(55, 1217-55: 1,2,1) 


С5, 1,1,1 
sinOcos8cosdsind 55, 2,1,2 
совбвіпфВ - sin0cos0sin?$ = = 258; „2.2.1 
cos20sinó$ 4С5, 11,1 


39 


E. PARTIALS OF W4 


аа - |А Гсов30віпфе(Е,0,ф)а8аф so the partials 


are given by: 


sin?6cos30cos?dsind 
2sin?0cos30cosósin?$ 
" 2віпбсовбсов3бсовфвіпф 
cosi r T sin?0cos30sin?$ 
25іпӨсоѕдсоѕзӨѕіп?ф 


сов2бсов3бвіпф 


2CS1,3,1,1 2051,3,3,1 CS1,5,1,1 *CS1,5,3,1 CS1,1,1,1 *CS1,1,3,1 
22СС,;,1,172СС з,3,1 7 СС15,1,17 СС 5,3,1 CC, 1, 1,1 * 06,3, 5,1) 
EE 55, 5,1,2 221,1,1,2 
24 2С5, ,з,1 CS1,5,3,1 CS1,1,3,1 
2(55;,5,2,17551,1,2,1) 


C51,5,1,13*2C59,,3,1,1* C91,4, 1,1 


ѕіпӨсоѕзӨсоѕзфѕіпф 
sinO0cos30sin?ó - 
cos0cos30sin$ 


5$51,4,1,2 991,2,1,2 
cos30sin$B = = 2(951,4,2,17991,2,2,1) 


2(C59, 4,1,1*C91,2,1,1) 
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Е. PARTIALS OF W5 


2 | . 
KESO aoo [ "f віп0еЕ,8,ф)40аф so the partials аге 
, , , 0 0 
glven by: 


sin*Ocos?o 
2віп?бсовфвіпф 
2sin?6cosb6cosd 
ѕіпӨА, = Е 
ѕіп°Өзіп?^ф 
2віп?бсовбвіпф 
cos20sin0 


5С; 1,2,1 
55, 1,1,2 
" 2(CC, 1,1,17 CC3,1,1,1) 
55, 1,2,1 
2(С5,|11-С5;114) 


SC, 1,0,0723,1,0,0 


sin?6cosb 25С, 11141 

š | ; 1. 
sin0B = ѕіп?Өѕіпф| = > 255, 1,1,1 
sin0cosO 5С 220,0 
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С. PARTIALS OF W6 


2u rn 
И = oc E f. jj sin20g(R,0,0)d0db so the partials are 
given by: 


sin?0sin20cos?$ 
2sin?Osin20cosdsind$ 
2sin0sin20cos0coso 
sin?*Osin20sin?$ 
2sinOsin20cosO0sind 
cos?0sin20 


sin20A, = 


25С, 2,1 9C1,4,2,1 
255 ;,1,27981 4,1,2 
(СС 1,17 СС,4,1,1525С6, 11,1) 
2551 2,2,1 991,4,2,1 ; 
2(С51 ,2,1,1—С51,4,1,112552,1,1,1) 
25С, 2,0,0* 95C1,4,0,0 


—- 
= 


E 
4 


sinO0sin20coso$ СС, 1,1,17 CC1,3,1,1 
sin20B = ѕ1п051п205іпф = = С 1,1, 17С а 


sin20cos0 SCA 3,0,0*9€4,1.0,0 
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H. PARTIALS OF W7 


2. . 
р боб - | [ віп300(Е,0Ө,ф)ідаф во «Пе раг“іа1в аге 
7 7 7 0 0 
given by: 


sin208sin30cos2ó 
2віп2бвіп3бсовфвіпф 
2sin0sin30cos0coso$ 

sin?Osin30sin?$ 
2sinO0sin30cos0sin$ 


sin30A, = 


cos?0sin30 


25С, ; ,175С11,2,175С15,2,1 
255, ; 1,2931 11,279515,1,2 
2(CC, 1,1,37 CC , 5,1,1) 
2551,3,2,1 991,1,2,1 $71,5,2,1 
2(CS1,1,1,2 CS1,5,1,1) 
SC1,5,0,0*2501,3,0,0*9C1,1,0,0. 


s 
4 


sin0sin30cosó$ СС, 2,1,1 CC1,4,1,1 
sinO0sin30sind| - = C51,2,1,1 C943,4,1,1 
sin30cos0 оо": об 


$1030 В 
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I. PARTIALS OF W8 


2X ^T 
Wi =n) Cy = f. iP sin*0cos$g(R,0,4)d0d$ so the partials are 


given by: 
ѕіп“Өсоз?ф 
2sin'6cos?bsind 
2sin?0cosOcos? 
ѕіп?ӨсоѕфА, = | ' N 
віп“Өзіп?2фсовф 
2віп?бсовбсовфвіпф) 
cos?0sin?0cos$ 


85C, 1,3,1 
16(95, 1,1,1^994,1,3,1) 
2(2 5C, 2,2,17 964, 4,2, 1) 
8(SC, 1,1,1^ 9C4,1,3,1) 
2551,2,1,2 SS1,4,1,2 
8(SC5 1,1,1^ 9C4,1,1,1) | 


= 1 
8 


L. 


25С, 12, 
55; 1,1,2 
2(CCi 1,1,37 CC3, 1,1, 1) 


sin*OÓcos?$ 
sin?0cos$B = |ѕіп?Өѕіпфсоѕф 
віп2бсовбсовф 


в" 
2 
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J. PARTIALS OF W9 


2х Л 
а 5С 7 аа J. js sin?6cos?$g(R,0,4)d0d$ so the partials are 


ѕіп°Өсоѕ“ф 
25іп?Өсоѕфѕ=іпф 
еси pul найи 
sin?*Osin?$cos?$ 
2sin*0cosOcos?$sin$ 
cos?0sin?0cos?$ 


4SCs 1,4,1 
555,1,1,4+2555,1,1,2 
8(СС, 1,з,1:2СС; 1,3,1%СС6 1,3,1) 
4(556 1,2,17556,1,4,1) 
8(CS1,1,1,1 CS1,1,3,1 2C95,1,1,1* 2C93,1,3,1* C9s,1,1,17 C98,1,5, 1) 
4(5С, 1,2,17566 1,2,1) 


— 


= + 
4 


sinf(0cos?$ 85C, 1,3,1 
віпЗбсов2фВ - |ѕіп*Өѕіпфсоѕ2ф| = =18(55,.1,1,175541,3,1) 
sin?0cosOcos?$ 25C, 2,2,1 9C1,4,2,1 
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К. PARTIALS OF W10 


W, 


p= SS aes ["Јзіпгөзіпфа, Ө, ф) аф во the partials are 


given by: 


віп“Өсов%фвіпф 

25іп“Өсовфвіп2ф 
а 2010 и 

sin'Osin'b 
2sin?0cosOsin?$ 
cos?0sin?0sinó$ 


8(55, 11,.17554 1,3,1) 
16(5С, 11175С, 1341) 
255, 1,2795314,1,2 


- i 
8 В 
455, :,172551 21 
8(55, 11,155, 1,1,1) 
сіпЗбсовфвіпф 55; 1,1,2 
sin?OsindB =| sin36sin?d | = z 255, 1,2,1 
sin*0cosO0sin$ 2(CS, 1,1,17 C95, 1,1, 1) 
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L. PARTIALS OF W11 


ТАИСА == fi” |] sin'esin26gR, 6, 6)d8d6 so the partials 


are given by: 


sin?O0cos?$sin24 
25іп?Өсоѕфѕіпфѕіп2ф 
обеси а "e Иа s 

ѕіп°?Өѕіп?фѕіп2ф 
2sin*0cosO0siné$sin26 
cos?0sin?0sin26 | 


SS5,1,1,4 *2555,1,1,2 
2(25C5,1,2,19Cs,1,1,4 SCs,1,1,2) 
a (CS11,1,3 + CS1,1,1,.1 20593,1,1,3 2093,1,1,2 *CS5,1,1,3 *CS5,1,1.1) 
2556 11,27956,1,1,4 
А(СС, 1,1,17СС61,1,1,3:2 СС; 1,1,152С С; 11,3%С61,1,17С66,1,1,3) 
4(55, 11,27956,1,1,2) 


= + 
4 


віп“Өсовфвіп2ф 4(SS, 1,1,3* 994,1,1, 1) 
sin'Osin268 - |sin*fOsinósin26| - 248644142561.) 
sin?0cosO0sin2$ 50,21 27951112 
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М. PARTIALS OF W12 


ZR PT a z 
Wie = ss sm J f. sin?Osin^$g(R,0,4$)d0db so the partials 
are given by: 


ѕіп°Өсоѕ2фѕіп2ф 
25іп?Өсоѕфѕіп?ф 
25іп“ӨсоѕӨсоѕфѕіп?ф 
ѕіп°Өѕіп“ф 
2sin*0cosOsin?$ 
cos?0sin?Osin?ó 


ѕіп?Өѕіп?фА, = 


4(55; 1:,1-556 1,4,1) 
255; 112-9556 11,4 


el (СС, 1,1,17С611,з,172С6; 11,1%2СС6; 1,3,1%С6,1,1,17С66,1,3.,1) 
4 455; 1,4,1 
8(CS, 1,3,172C95,1,3,1* C95,1,5,1) 
| 4(55; 1,2,17556,1,2,1) 
ѕіп*Өсоѕфѕіп?ф 8(SC, 1,1,17 9C4,1,3,1) 
sin?Osin?$B - sin*Osin*$ = = 855, 1,3,1 
sin'Ocos0sin?$ 255, 2,2,1 7 991,4,2,1 
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APPENDIX C 
A. PARTIALS OF U1,U2 AND U3 


This step simply involves taking derivatives of 


[Ur u, u,]7 E 








T 
Wa Wio We | 
Ws М 2W, 


with respect to each parameter. The results are: 











ди, w 
Ow. м2 
Ca Т 
Ow, W; 
ди, mio 
Ow. we 
OD ТТ 
ди, W, 
д», 22 
ou 1 
Әу 2 


and all other partials = 0. 
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В. PARTIALS OF U4,..,U9 WRT MU 














E 
Qu, O W| 

9 2 W, 

P zs 
9и, = ag g12 
W Ms | 13 

G 
pan 
Qu W, 
a = 6 gi? 
H ES 13 
G 
12 
сй а EU UE 
EE 
P piis 23 
G 
16 
див С 2 
9 2 w. 
P 21225 
B 
Qu и 

9 
5 qum 6 923 

P 2” 33 

G 
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С. PARTIALS OF U4,..,U9 WRT SIGMA 
To solve for these partials we use the table of partials 


given in the previous section and use standard calculus rules: 
1. Partials of U4 






































our колос Ки: Ю.и шо Ки, — wo) _ Изв 
до,» ОО Г mw 271775) ву 2 ws 
du, E om 1 (g Ки, -w,) MHM3Wg|- Ки: 829% 
00,, 4 gq?| "| 8w 2 ws GM MEN Ww. 
gu 6224-1, R(W,-W,) _ вм. \_ Rw, MWg 
DOM q2 (| 8w, 2 Ws 9. 2и; 2м; 
ди, = -2o23u,-— |o, rum 
4 3 
до.» 2 Ws 
D Ки: а к” R(W,-W;) M3Wg 
mM Iw NEP вк 2 ws 
Qu, 








E» 


2. Partials of U5 





























ди Rw,,-24.w. 
= = -20u i lon ee) 
023 q Ws 
_ Ки |} Қ(м;-%) о 
оз ЕЕ КӨЗ: теи сти 
Ws Ws Ws 
du, _ E Al a Ки, 0 d RW,  MiWig 
дозз а?| A м; Ws 22 2%5 Ws 





3. Partials of U6 
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Ou, 
































= ТИЕ [<> S Es W Ки. -и.) _ Иа 
don. ° < E) nos 4 ws 2 ws 
ди, к || ae LAA E: R(W,-W,)) р, м, 
DOM 6 g? A AU IZ | aw, 2%; 
Qu, = -2673yu,-+ a RWi-W;) 4 
005, а? 2%; "т 
is еее СЕКЕ "s |-o d ach Ws) _ Я] 
13 
4 w; 2 w; 2 w; 
du, TONES M R(W,-W,) H2% n R(W,7W;) bi We 
005; t c МЧ» 2w,} 74 ам, 2%; 
4. Partials of U7 
ди, E ашна R(W,-W,) M3Wio Ж RW Hio 
95.1 7 g?| "| 8w 2w. 33 2. 2w, 
du, — ао... „ын Ки Mit) К(Мҙ-М) _ Ше 
00,, 7 q2 U aw. 2», 913 Sw. 2м 
du, Sy ede ele Rw Hio 
Оо,» Ge Ws 


КУ 1 Кизи.) -W 4) _ 


8w- 


ed W за ые 


4 w; 
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Haw 


4 


2 w. 





ди, - -2o02u,- 1. i. | Wi) Ha а). N Ri | 
g? 


Ws 2W 4 W, 2 Ws 











du, = -0u NM RW,, И.и ES Ки; H2Wio 
Qo 7 g?| 74. 4w 2 w, Ee E 2 Ww, 
33 q 5 5 5 e i 


5. Partials of U8 


Qus = -glly E (522%. a [35 )- | 
8 = 23 "SP 33| eS 


00i, 4 Ws 2 W; 4 ws 2 We 





Ou, = -20124,--- a 6, am R(W1-Ww,) -%;)- tt) ono R(w,*tw.) s 


do, 4 ws 2 W; 4 ws 2 w; 

















du, - -209u,- 1. i K(W,7W,) M2W; 
00,4 = |0 2 и. Ws 
| 2 wa) -W,) MiWe R(W,tWs) |; 
4w 2 Ww, 1 ди. 2%. 
диз _ — 9221) 
005, 
Ou, _ 20231) 5“. W,-W,) 51% 
00,5 ІІ AW. “ 4 2, 


5а 




















6. Partials of U9 
Qu, — ЧОП ЖЕНЕ ЕЕЕ 
do Чә 2| 23 8 
Ti q : 
ди R(w =) 
2 = -202u,-— i 7 
до, g? 
Е R(W3-W,) HW 
Orc 
5 
ди W, -W 
J а= -20%3u,-— 0,, RW Wa) 
90; аг 8и. 
ди. = -g?? Mr EL Е W. W, 
do > 13 a 
22 q? 
du, E E W.W.) 
до. g? 8 w. 
ди» 
00,, 
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0; Киги) В: -0,, Киз -и,) 
D 4 W, 4 Ww 


_ №21 
2 Ws 


2 Ws 


ioc -g KW] + Ws) N 
4 W, 22 Sw 
5 5 


E 
2 Ws 


H.W, 
4 Ws 





В.и |_ R(W,-W,) 
О: <. 
4 Ws 8 w; 
13% 
4 W. 


к: B, Wç 
4 Ws 


x [Хеге 


4 W; 8 w; 


-0 R(w,-W3) 
a т 12 8 Ws 


D. PARTIALS OF U4,..,U9 WRT R 














A i ТАНА ТЫЛУЫ) 

=: 2 qu, 07; 401 Qus -w) 

95 Ë ар (2077 Ha *4022w,, +o22(w, -w.)) 
A | au o 0 79) * oH s cv) o?) 


E. PARTIALS OF U4,..,U9 WRT W 


1. Partials of U4 


ди. on Roi? 
Ow, 8 w, 

Qu, E Rg?? 
Ow, 8 w, 
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2. 














ди 
4 = EU 
ди. W5 
ди, l 
Pn 1 2 
Ow, | (o uto uto y) 
Qu, Роіі 
Ow, 2 W, 
ди, Rol 
OW,, 4 W. 
Partials of U5 
Ow, 4 W. 
OW, ауу 
Qus 1 
- — us 
OW. Ws 
Qus 
ВИ 11 12 13 
OW, < mm В рз) 
ди» _ Ра11 
Әу, 2 W. 
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3. Partials of U6 


ди; _ Rott 

















Ow, 4 W; 
ди; _ Roll 
9%, 4 Ws 
Qu, Ro12 
Ow, 4 Ws 
9% _ Rol2 
ду, 4w 
13 
ды _ mu + 
Ow, Ws 4 W. 
Qu, 1 
2D 11 +G 12 +013 
OW, 25; EON 


Qu, коз 
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4. Partials of U7 





























Ou, а Rg?? 
ду, 8 w; 
ди, _ Rg?? 
Ow, 8 w; 
ди 
EA = l U, 
Ow, Ws 
du, 1 
Кее оь) 
ди, Е Roi? 
ду, 4w; 
dU, _ Ro? 
ду’, > 2 w; 
5. Partials of U8 

ди, _ Ror 
Ow, 4 w. 
ди; Rot 
Ow, ау; 

du, W Ro??? 

Ow, AW. 
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Ow, 4 Ws 
ди; _ 1 Ro’? 
ди, 
c сі? +022 +923 
Ow, Бы. E a s 


6. Partials of 09 


ди, _ Rolt? 














Ow, 8 Ws 
Ou, _ Ro? 
Ow, 8. 
ди, _ Ro23 
Ow, 8 w; 
Ou, _ _ Ro23 
Ow, 8 Ws 
OU LE Ro 
OW. Ws i 8 Ws 
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1 
B (op, to toy) 





Ow, 4 w; 
ди, Rg?? 
Ow, 8 w, 
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* XM 0X 0X f * X 


APPENDIX D 


PROGRAM SEPCMP 


PROGRAMMER : LT ARTHUR F. BROCK 
DATE: 12 MAY 1991 
LAST MODIFIED: 29 AUGUST 1991 
PROGRAM DETERMINES FIRST AND SECOND ORDER APPROXIMATIONS TO 


VARIANCE OF SEP GIVEN VALUES OF MU AND SIGMA (VARIANCE - COVARIANCE 
MATRIX) AS INPUTS. OTHER INPUTS ARE PARAMETERS FOR USE IN THE 


IMSL SUBROUTINE QTWODQ USED TO PERFORM GAUSSIAN QUADRATURE 


EVALUATION OF DOUBLE INTEGRALS. 


20 


10 


50 
51 
52 


53 
55 


INCLUDE !COM DEF! 
INCLUDE 'WVEC DEF! 
INCLUDE 'PMUCOM DEF! 
INTEGER 1,011,072, SAMPSZ,K 
REAL DT1(10,9),DT2(22,10),DT3(9,22),SEPMC(9,9),PMUC(3,3),PSIG(6,6) 
, SEPMUC(3),SEPSIG(6),CI(2),MPT(3,3) 
CALL INIT(SAMPSZ,MPT) 
WRITE(16,51) R,NUMTR 
WRITE(16,52)(MU(1),121,3) 
DO 20 121,3 
WRITE(16,53) (MSIG(I, J), J21,3) 
CONTINUE 
OT 1=0 
OT2=0 
DO 10 1=1,SAMPSZ 
CALL INIT2(MPT) 
CALL WVALS 
CALL DT1VAL(DT1, SEPMU, SEPSIG) 
CALL DT2VAL(DT2) 
CALL DT3VAL(DT3) 
CALL SEPEST(DT1,DT2,DT2, SEPM) 
CALL PMUS(MSIG,PMU,PSIG) 
CALL BNDEST(SEPM,PMU, PSIG, SEPMU, SEPSIG,CI) 
IFCSEPSET .LT. (R-CI(1)) .OR. SEPSET .GT. (R*CI(1))) THEN 
OT1=0T 1+1 
IFC(SEPSET .LT.(R-CI(2)).0R. SEPSET .GT.(R*CI(2))) THEN 
OT2=0T2+1 
ENDIF 
ENDIF 
PRINT*,'I = ',1,REALCI-OT1)/REAL(1),REALCI-OT2)/REAL(1) 
WRITEC16,50) R-CI(1),R+CI(1),0T1,R-CIC2),R+CI(2),012 
CONT I NUE 
WRITE(16,55) SAMPSZ,REAL(OT1)/REAL(SAMPSZ) ,REAL(OT2) /REAL(SAMPSZ) 
FORMAT(1X,20' —(!',F7.2,',',F7.2,')',2X,13,2X)) 
FORMAT(1X,'TEST FOR SEP = ',F7.2,/,1X, 'SAMPLES PER TRIAL = ',12) 
FORMAT(1X, ! MU ',/,3CF7.2,',"),/,/, 1X, 


E ; SIGMA') 


FORMAT(1X,3(F9.3,3X)) 


FORMAT(1X,/,'SAMPLE SIZE: ',14,3X,2(F5.3,8X)) 
STOP 
END 
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SUBROUTINE INITCSAMPSZ,MPT) 
INITIALIZE FILES, READ INPUT DATA AND INITIALIZE VARIABLES. 


INCLUDE 'COM DEF! 

INCLUDE 'PMUCOM DEF! 

INTEGER N,SAMPSZ,RK 

REAL MPT(3,3),TOL,NUMSIG(3,3) 
EXTERNAL RNSET,DCHFAC 
OPEN(15,FILE="!/MUSIG DATA Al!) 
OPEN( 16, FILE="/OUTDT DATA A1') 


READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15, *) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 
READ(15,*) 


MU( 1) 
MU(2) 
MU(3) 
SIG(1) 
SIG(2) 
SIG(3) 
SIG(&) 
SIG(5) 
516(6) 
ERRABS 
ERRREL 
IRULE 
NUMTR 
SAMPS2 


READ(15,*) ISEED 

READ(15,*) NUMSIM 

CALL RNSET(ISEED* INT(MU(1)*SIG(1))) 
DET=SIG(1)*(SIG(4)*SIG(6)-(SIG(5)**2))-SIG(2)*(SIG(2)*SIG(6) 
C -SIG(3)*SIG(5))*SIG(3)*(SIG(2)*SIG(5)-SIG(3)*SIG(4)) 
IF (DET .LE. 0.0000000001) THEN 

WRITE(16,50) DET 
STOP 

ENDIF 
SIGINV(1)=(SIG(4)*SIG(6)-(SIG(5)**2))/DET 
SIGINV(2)=(SIG(3)*SIG(5)-SIG(2)*SIG(6))/DET 
SIGINV(3)=(SIG(2)*SIG(5)-SIG(3)*SIG(4))/DET 
SIGINV(4)=(SIG(1)*SIG(6)-(SIG(3)**2))/DET 
SIGINV(5)=(SIG(3)*SIG(2)-SIG(1)*SIG(5))/DET 
SIGINV(6)=(SIG(1)*SIG(4)-(SIG(2)**2))/DET 
DENOM=(SORT(2.0* PI )**3)*SORT(DET) 
М516(1,19-516(1) 

MSIG(1,2)=SIG(2) 

MSIG(1,3)=SIG(3) 

MSIG(2,1)=S1G(2) 

М516(2,2)=516(4) 

М516(2,3)=516(5) 

М516(3,1)=516(3) 

М$16(3,2)=$16(5) 

М$16(3,3)=$16(6) 

CALL FINDSEP 

MSIG(1,1)=SIG(1)/NUMSIM 
MSIG(1,2)=SIG(2)/NUMSIM 
MSIG(1,3)=SIG(3)/NUMSIM 
MSIG(2,1)=SIG(2)/NUMSIM 
MSIG(2,2)=SIG(&)/NUMSIM 
MSIG(2,3)=SIG(5)/NUMSIM 
MSIG(3,1)=SIG(3)/NUMSIM 
MSIG(3,2)=SIG(5)/NUMSIM 
MSIG(3,3)=SIG(6)/NUMSIM 

TOL=100.0*DMACH(4) 

CALL DCHFAC(3,MSIG,3,TOL,RK,MPT,3) 

CALL SETDAT 

CALL MAKES 

FORMAT(1X,' ILLEGAL INPUT MATRIX, DETERMINANT = ',F8.5,' PROGRAM', 
C ' TERMINATED. ') 

RETURN 

END 
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20 


15 


30 


20 
10 


SUBROUTINE INIT2(MPT) 
INITIALIZE VALUES USED EACH TIME THROUGH LOOP. 


REAL MPT(3,3),TMU(3),TSIG(Ó) 
INCLUDE 'COM DEF' 
INCLUDE 'PMUCOM DEF' 
CONTINUE 
$16(1) 
SIG(2) 
SIG(3) 
SIG(4) 
SIG(5) 
516(6) 
MU(1)20. 
М0(2)-0.0 
М0(3)-0.0 
DO 15 1=1,INT(NUMSIM) 
CALL TRIALSC(NUMTR,MPT, TMU, TSIG) 
SIG(1)=SIG(1)+TSIG(1) 
SIG(2)=SIG(2)+TSIG(2) 
SIG(3)=SIG(3)+TSIG(3) 
SIG(4)=SIG(4)+TSIG(4) 
$16(5)=$16(5)+1$16(5) 
SIG(6)=SIG(6)+TSIG(6) 
MU(1)2MUC1)* TMUC1) 
MUC(2)=MU(2)+TMU(2) 
MU(3)=MUC3)+TMUC3) 
CONTINUE 
М516(1,19-516(1) 
MSIG(1,2)=SIG(2) 
MSIG(1,3)=S1G(3) 
MSIG(2,1)=SIG(2) 
MSIG(2,2)=SIG(4) 
MSI G(2,3)=SIG(5) 
MSIG(3,1)=SIG(3) 
MSIG(3,2)=SIG(5) 
MSIG(3,3)=SIG(6) 
DET=SIG(1)*(SIG(4)*SIG(6)-(SIG(S)**2))-SIG(2)*(SIG(2)*SIG(6) 
-$16(3)*$16(5))+$16(3)*($16(2)*$16(5)-$16(3)*$16(4)) 
IF (DET .GT. 0.00000000001) GOTO 30 
GOTO 20 
CONTINUE 
SIGINV(1)=(SIG(4)*SIG(6)-(SIG(5)**2))/DET 
SIGINV(2)=(SIG(3)*SIG(5)-SIG(2)*SIG(6))/DET 
SIGINV(3)=(SIG(2)*SIG(5)-SIG(3)*SIG(4))/DET 
SIGINV(4)2(SIG(1)*SIG(Ó) - (SIG(3)**2) )/DET 
SIGINV(5)2(SIG(3)*SIG(2)-SIG(1)*SIG(5))/DET 
SIGINV(6)=(SIG(1)*SIG(4)-(SIG(2)**2))/DET 
DENOM-(SQRT(2.0*PI)**3)*SQRT(DET) 
CALL FINDSEP 
RETURN 
END 


оооооо 


О • 
EO 000.00 


SUBROUTINE MXTMLT(M1,M2,ANS,1,J,K) 


PERFORMS MATRIX MULTIPLICATION M1(1,J)*M2(J,K) TO PRODUCE 
OUTPUT MATRIX ANS(CI,K) 


INTEGER L,ROW,COL,1,J,K 
REAL M1C1,J),M2CJ,K),ANSC1,K) 
DO 10 ROW=1,1 
DO 20 COL=1,K 
ANS(ROW,COL)=0.0 
CONTINUE 
CONTINUE 
DO 30 ROW=1,1 
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DO 40 COL=1,K 
DO 50 L=1,J 
ANS(ROW, COL )=ANS(ROW, COL)+M1(ROW,L)*M2(L, COL) 
CONT I NUE 
CONT I NUE 
CONT I NUE 
RETURN 
END 


SUBROUTINE SEPEST(DT1,DT2,DT3, SEPM) 
PRODUCES THE 9 BY 9 MATRIX OF PARTIALS OF SEP 


INCLUDE 'COM DEF! 

REAL DT1(10,9),DT2(22,10),0T3(9,22), SEPM(9,9),TEMP(22,9) 
CALL MXTMLT(DT2,DT1, TEMP, 22, 10,9) 

CALL MXTMLT(DT3, TEMP, SEPM,9, 22,9) 

RETURN 

END 


SUBROUTINE SETDAT 


SETS THE INPUT VALUES OF MU, SIG AND SEP TO BE HELD CONSTANT 
THROUGHOUT THE LOOPS OF THE PROGRAM. 


INCLUDE 'COM DEF! 
INCLUDE 'PMUCOM DEF! 
INTEGER I,J 
0010 1-1,3 

MUSET(I )=MUCI) 
CONTINUE 
SEPSET=R 
RETURN 
END 


SUBROUTINE TRIALS(N,MPT,TMU,TSIG) 


DRAWS N RANDOM SAMPLES FROM NORMAL(MU,SIGMA) AND DETERMINES 
ESTIMATES FOR MU AND SIGMA BASED ON THESE SAMPLES. 


INCLUDE 'COM DEF' 

INCLUDE 'PMUCOM DEF! 

INTEGER I,N,J 

REAL MPT(3,3),RVAR(50,3),RSQ(50,3),RMLT(50,3),RMN(6),OP(3) 

,NUMFAC, TSIG(6), TMUC(3) 

EXTERNAL DRNMVN, RNSET,DMACH, DRNUNF 

00 10 1-1,М 
CALL DRNMVN(1,3,MPT,2,0P, 1) 
RVARCI , 1)20P( 1) *MUSET(1)/NUMSIM 
RVAR(1, 2)=OP(2)+MUSET(2)/NUMSIM 
RVARC1,3)=OP(3)+MUSET(3)/NUMSIM 
RSQ(I, 1)-RVAR(CI, 1)**2 
RSQ(I,2)-RVAR(I,2)**2 
RSO(1,3)=RVAR(1,3)**2 
RMLT(I, 1)-RVARCI , 1)*RVARCI,2) 
RMLT(I,2)-RVARCI , 1)*RVAR(CI,2) 
RMLT(1,3)=RVAR(I,2)*RVARCI,3) 


10 CONTINUE 


TMUC1)=FMNCRVAR,1,WN) 
TMU(2)=FMN(RVAR,2,N) 
TMU(3)=FMN(RVAR,3,N) 
RMN(1)=FMNC(RSQ,1,N) 
RMN(2)=FMN(RSQ,2,N) 
RMN(3)=FMN(RSQ,3,N) 
RMN(4)=FMNC(RMLT, 1,N) 
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RMN(5 )=FMNCRMLT,2,N) 
RMN(6)=FMNCRMLT,3,N) 

NUMFAC=1.0 

TSI1GC1)z CRMNC1)- TMUC1)**2) *NUMFAC 
TSIG(2)=(RMN(4)-TMUC1)*TMUC2))*NUMFAC 
TSIG(3)=(RMNC5)-TMUC1)*TMUC(3))*NUMFAC 
TS$1GCA)  CRMN C2) - THUC2) **2) *NUMFAC 
TSIG(5)=(RMN(6)-TMU(2)*TMU(3))*NUMFAC 
TS$1GC 6) CRMN (3) - THUC3) **2) *NUMFAC 
RETURN 

END 


REAL FUNCTION FMN(VEC,COL , N) 
DETERMINES THE MEAN OF N ITEMS IN COLUMN COL OF VECTOR VEC 


INTEGER I,N,COL 

REAL VEC(50,3),TOT 

TOT=0.0 

00 10 1-1,М 
TOT=TOT+VEC(1, COL) 

CONT I NUE 

FMN=TOT/REAL(N) 

RETURN 

END 


REAL FUNCTION G(X) 
REAL X 

G=0.0 

RETURN 

END 


REAL FUNCTION H(X) 
REAL X 

Н=6.2831854 

RETURN 

END 


SUBROUTINE WVALS 


DETERMINES VALUES OF W, CC, CS, SC AND SS FUNCTIONS USING THE 
IMSL SUBROUTINE QTWOOQ. 


REAL G,H,H1,H2,ERREST 

INCLUDE 'COM DEF! 

INCLUDE 'WVEC DEF! 

EXTERNAL G,H,DTWODQ, CC1111,CC1211,CC2111,CS1111,CS1211,CS2111, 
5С1100,5С1200,5С1300,5С2111,5С3121,552111,553112,553121, 
СС1113,СС1131,СС1211,СС1331,СС1411,СС1511,СС1531,СС2111, 
СС3113,СС3131,СС5111,СС5113,СС5131,С51113,С51131,С51211, 
С51331,С51411,С51511,С51531,С52111,С53113,С53131,С55111, 
С55113,С55131,5С1121,5С1212,5С1221,5С1321,5С1400,5С1412, 
5С1421,5С1500,5С1521,5С3100,5С4111,5С4113,5С4131,5С5112, 
5С5114,5С5121,5С5141,551112,551121,551212,551221,551312, 
551321,551412,551421,551512,551521,554111,554113,554131, 
555112,555114,555121,555141 

В=РІ 

CALL DTWODQ (CC1111,A,B,G,H,ERRABS, ERRREL, IRULE,W(1),ERREST) 

CALL DTWODQ (CC1311,A,B,G,H, ERRABS, ERRREL, IRULE,W(2), ERREST) 
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(SC5141, 
ERRABS, ERRREL , IRULE , SS( 1) , ERREST) 
ERRABS, ERRREL , IRULE, SS(2) , ERREST) 
ERRABS, ERRREL, IRULE, SS(3), ERREST) 
ERRABS, ERRREL , IRULE , SS(4) , ERREST) 
ERRABS, ERRREL, IRULE,SS(5),ERREST ) 


(CS1111,A,8,G,H, ERRABS, ERRREL, IRULE,W(3), ERREST) 
(CS1311,A,B,G, H, ERRABS, ERRREL , IRULE , W(4) , ERREST) 
(SC1100,A,B,G, H, ERRABS, ERRREL, IRULE , W( 5) ,ERREST) 
(SC1200,A,B,G, H, ERRABS, ERRREL, IRULE , W( 6) , ERREST) 
(SC1300,4,B,G,H,ERRABS, ERRREL, IRULE,W(7), ERREST) 
(SC2111,A,B,G, H, ERRABS, ERRREL, IRULE, W(8) , ERREST) 
(SC3121,A,B,G, H, ERRABS, ERRREL, IRULE , W(9) , ERREST) 
(SS2111,A,B,G, H, ERRABS, ERRREL , IRULE , W( 10) , ERREST) 
(SS3112,A,B,G, H,ERRABS, ERRREL , IRULE, W( 11) , ERREST) 
(SS3121,A,B,G, H, ERRABS, ERRREL , IRULE , W( 12) ,ERREST) 
(CC1113,A,B,G, H, ERRABS, ERRREL , IRULE , CC( 1) , ERREST) 
(CC1131,A, B, G, H, ERRABS, ERRREL , IRULE , CC( 2) , ERREST) 
(CC1211,A,B,G, H, ERRABS, ERRREL,, IRULE , CC( 3) , ERREST) 
(CC1331,A,8,G,H, ERRABS, ERRREL, IRULE, CC( 5), ERREST) 
(CC14511,4,B, G, H, ERRABS,, ERRREL , IRULE , CC( 6) , ERREST) 
(CC1511,A,B,G,H, ERRABS, ERRREL, IRULE, CC(8), ERREST ) 
(CC1531,A,B, G, H, ERRABS, ERRREL , IRULE, CC(9) , ERREST) 
(CC2111,A, B, G, H, ERRABS, ERRREL , IRULE , CC( 10) ,ERREST) 
(CC3111,A,B, G, H, ERRABS, ERRREL,, IRULE , CC( 17), ERREST) 
(CC3113,A,B,G, H, ERRABS, ERRREL, IRULE , CC( 11) ,ERREST) 
(CC3131,A,B, G, H, ERRABS, ERRREL , IRULE , CC( 12) , ERREST) 
(CC5111,4,B,G, H,ERRABS, ERRREL , IRULE, CC( 14) , ERREST) 
(CC5113,A,B,G,H, ERRABS, ERRREL, IRULE,CC( 15), ERREST) 
(CC5131,A.B,G,H, ERRABS ERRREL , IRULE , CC( 16) , ERREST) 
(CS1113,A,B,G, H, ERRABS, ERRREL , IRULE, CS(1) , ERREST) 
(CS1131,A,B,G, H, ERRABS, ERRREL , IRULE, CS(2) , ERREST) 
(CS1211,A,8,G,H, ERRABS, ERRREL, IRULE, CS(3), ERREST) 
(CS1331,A,B,G,H, ERRABS, ERRREL, IRULE,CS(5), ERREST) 
(CS1411,A,B,G, H, ERRABS, ERRREL, IRULE, CS(6), ERREST) 
(CS1511,A,B,G,H, ERRABS, ERRREL, IRULE, CS(8), ERREST) 
(CS1531,A,B,G, H, ERRABS, ERRREL , IRULE , CS(9) , ERREST) 
(CS2111,A,B, G, H, ERRABS, ERRREL , IRULE , CS(10) , ERREST) 
(CS3111,A,B,G, H, ERRABS, ERRREL , IRULE , CS(17) , ERREST) 
(CS3113,4,B,G, H, ERRABS, ERRREL , IRULE, CS( 18) , ERREST) 
(CS3131,A,B,G, H,ERRABS,ERRREL , IRULE,CS( 11) , ERREST) 
(CS5111,A,B,G, H,ERRABS, ERRREL,, IRULE, CS( 13) , ERREST) 
(CS5113,A,8,G,H, ERRABS, ERRREL, IRULE, CS( 14) , ERREST) 
(CS5131,A,B,G, H, ERRABS, ERRREL, IRULE , CS( 16) , ERREST) 
(SC1121,A,8,G, H, ERRABS, ERRREL , IRULE , SC( 1) , ERREST) 
(SC1212,A,B,G, H, ERRABS, ERRREL , IRULE, SC( 18) , ERREST) 
(SC1221,A,B,G, H, ERRABS, ERRREL , IRULE , SC(2) , ERREST) 
(SC1321,A,B,G,H, ERRABS, ERRREL , IRULE , SC(3) , ERREST) 
(SC1400,A,B,G, H, ERRABS, ERRREL, IRULE , SC(4) , ERREST) 
(SC1412,A,B,G, H, ERRABS, ERRREL, IRULE , SC( 19) , ERREST) 
(SC1421,A,B,G, H, ERRABS, ERRREL, IRULE , SC(5) , ERREST) 
(SC1500,4,B,G, H, ERRABS, ERRREL , IRULE , SC(6) , ERREST) 
(SC1521,A,B,G, H, ERRABS, ERRREL , IRULE , SC(7) , ERREST) 
(SC3100,4,B,G, H, ERRABS, ERRREL , IRULE , SC(9) , ERREST) 
(SC4111,A,B,G, H, ERRABS, ERRREL , IRULE , SC( 10), ERREST) 
(SC4113,A,B,G,H, ERRABS, ERRREL, IRULE, SC( 11), ERREST) 
(SC4131,A,B,G, H, ERRABS, ERRREL , IRULE, SC( 12) , ERREST) 
(SC5112,A,B,G, H, ERRABS, ERRREL, IRULE, SC( 13) , ERREST) 
(SC5114,A,B,G, H, ERRABS, ERRREL, IRULE , SC( 14) , ERREST) 
(SC5121,A,B, G, H, ERRABS, ERRREL , IRULE , SC( 15) , ERREST) 

A,B,G 

A,B,G 

A,B,G 

A,B,G 

A,B,G 

A,B,G 

A,B,G 

‚А,В,С 

А,В,С 

„А,В,С 

‚А,В,С 

А,В,С 

А,В,С 

А,В,С 


, 
, 
, 
, 
, 
, 
, 
, 
, ERRABS, ERRREL, IRULE, SC( 16) , ERREST) 
, 
, 
, 
, 
, 


(SS1321,A,B,G,H,ERRABS,ERRREL,IRULE,SS(6),ERREST) 
(SS1412,A,B,G, H, ERRABS, ERRREL, IRULE ,SS(7), ERREST) 
(SS1421,A,B,G, H, ERRABS, ERRREL , IRULE, SS(8) , ERREST) 
(SS1512,A,B,G, H, ERRABS, ERRREL,, IRULE , SS(9) , ERREST) 
(SS1521,A,B,G, H, ERRABS, ERRREL , IRULE, SS( 10) , ERREST) 
(SS4111,A,B,G, H, ERRABS, ERRREL, IRULE, SS( 12) , ERREST) 
(SS4113,A,B,G, H, ERRABS, ERRREL , IRULE, SS( 13) , ERREST) 


(SS4131, ERRABS , ERRREL , IRULE, SS( 14) , ERREST) 
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CALL DTWODQ (SS5112,A,B,G, H, ERRABS, ERRREL, IRULE, SS( 15) ,ERREST) 
CALL DTWODQ (SS5114,A,B, G, H, ERRABS, ERRREL, IRULE, SS( 16) , ERREST) 
CALL DTWODQ (SS5121,A,B,G, H, ERRABS , ERRREL , IRULE , SS(17) ,ERREST) 
CALL DTWODQ (SS5141,A,B,G, H, ERRABS, ERRREL , IRULE, SS(18) , ERREST) 
RETURN 

END 


SUBROUTINE FINDSEP 


BASED ON PL-1 PROGRAM RAP, PRODUCES SEP BASED ON MU AND SIGMA 
FINDSEP DIAGONALIZES SIGMA BASED ON EIGENVALUES AND THEN CALLS 
SUBROUTINE ITERATION, WHICH DETERMINES SEP 


INCLUDE 'COM DEF' 
REAL EIGVLS(3),EIGMAT(3,3),01(3,3) , NMEAN(3) , NEMSIG(3, 3) 
C — ,TEMP(3,3),EIGNML(3,3) 
INTEGER I,J 
EXTERNAL DEVCSF 
CALL DEVCSF(3,MSIG,3,EIGVLS,EIGMAT, 3) 
CALL ORTHO(EIGMAT, EIGNML) 
DO 11 1=1,3 
DO 21 J=1,3 
Q1(I,J)sEIGNML(J, I) 
CONTINUE 
CONTINUE 
CALL MXTMLT(Q1,MU,NMEAN,3,3,1) 
CALL MXTMLT(MSIG,EIGNML, TEMP,3,3,3) 
CALL MXTMLT(Q1,TEMP,NEWSIG,3,3,3) 
CALL ITERATION(NMEAN, NEWSIG) 
RETURN 
END 


SUBROUTINE ORTHO (A1,A2) 
ORTHONORMALIZES MATRIX OF EIGENVECTORS OF SIGMA 


REAL A1(3,3),A2(3,3),L, TEMP, TEMP2 
INTEGER I 
L=SORT((A1(1,1)**2)+(A1(2,1)**2)+(A1(3,1)**2)) 
DO 10 1=1,3 

A2(1,1)=A1(1,1)/L 
CONTINUE 
ТЕМР-А2(1,12%А1(1,2)%А2(2,1)%А1(2,2)%А2(3,1)%А1(3,2) 
00 20 1-1,3 

А2(1,2)-А1(1,2)-ТЕМР%А2(1,1) 
CONT INUE 
L=SORT(A2(1,2)**2+A2(2,2)**2+A2(3,2)**2) 
DO 30 1=1,3 

A2(1,2)-A2(1,2)/L 


50 CONTINUE 


40 


TEMP=A2(1,1)*A1(1,3)+A2(2,1)*A1(2,3)+A2(3,1)*A1(3,3) 
TEMP2=A2(1,2)*A1(1,3)+A2(2,2)*A1(2,3)+A2(3,2)*A1(3,3) 
DO 40 1=1,3 
А2(1,3)-А1(1,3)-ТЕМР%А2(1,1)-ТЕМР2%А2(1,2) 
CONTINUE 
L=SORT(A2(1,3)**2+A2(2,3)**2+A2(3,3)**2) 
DO 50 1=1,3 
A2(1,3)=A2(1,3)/L 


50 CONTINUE 


RETURN 
END 
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SUBROUTINE ITERATION(MEAN, SIGMA) 


GENERATES SEP BASED ON MU AND SIGMA FOR DIAGONALIZED SIGMA 


INCLUDE 'COM DEF' 
INCLUDE 'SEPCOM DEF' 
REAL TOL,MEAN(3),SIGMA(3,3),MOM1, TRACE(3),SQSIG(3,3),CUSIG(3,3) 


С , TEMP(3),C2,D0F,RLOW, TPLOW, RHIGH, TPHIGH, CRATIO, INTEGRAL 
С , MOM2 , MOM3 
TOL=0.001 


СХХ-1.0/516МА(1,1) 
CYY=1.0/SIGMA(2, 2) 
CZZ=1.0/SIGMA(3, 3) 
CSANT=QUP I TW*SQRT(CXX*CYY*CZZ) 
C=(CXX* (MEAN(1)**2)+CYY*(MEAN(2)**2)+CZZ*(MEAN(3)**2))/2.0 
CALL MXTMLT(MEAN,MEAN,MOM1,1,3,1) 
CALL MXTMLT(SIGMA, SIGMA, SQS1G, 3, 3,3) 
CALL MXTMLT(SIGMA,MEAN,TEMP,3,2, 1) 
CALL MXTMLT (MEAN, TEMP,MOM2, 1,3, 1) 
CALL MXTMLT(SIGMA,SQSIG,CUSIG,3,3,3) 
CALL MXTMLT(SQSIG,MEAN, TEMP,3,3, 1) 
CALL MXTMLT (MEAN, TEMP,MOM3,1,3,1) 
00 10 1-1,3 
ТРАСЕ(1)-0.0 
MN(I)=MEAN(I) 
CONT INUE 
00 20 1-1,3 
ТРАСЕ(1)-ТВАСЕ(1)%516МА(1,1) 
TRACE(2)=TRACE(2)+SQSIG(I,1) 
TRACE(3)=TRACE(3)+CUSIG(I,1) 
CONTINUE 
MOM 1=MOM1+TRACE (1) 
MOM2=2.0* (TRACE (2)+2.0*MOM2) 
MOM3=8 .0* ( TRACE (3)+3 .0*MOM3) 
BETA=(MOM3**2)/(MOM2**3) 
DOF=8.0/BETA 
C2-DOF*(1.0-(2.0/(9.0*DOF)))**3 
RAD IUS=SORT (ABS(SORT(ABS(MOM2/(2.0*DOF)))*(C2-DOF )+MOM1)) 
RLOW=0.95*RADIUS 
TPLOW=EVALT(RLOW) 
CONTINUE 
IF(TPLOW .GT. 0.5) THEN 
RLOW=0.9*RLOW 
TPLOW=EVALT (RLOW) 
GOTO 40 
ENDIF 
RHIGH=1. 105*RLOW 
TPHIGH=EVALT(RHIGH) 
CONTINUE 
IF (TPHIGH .LT. 0.5) THEN 
RHIGH=1.1*RHIGH 
TPHIGH=EVALT(RHIGH) 
GOTO 50 
ENDIF 
CRAT10=(0.5-TPLOW)/(TPHIGH-TPLOW) 
RAD IUS=RLOW+(RHIGH-RLOW)*CRATIO 
ICNT=0 
CONTINUE 
ICNT=ICNT+1 
IF((RHIGH-RLOW)/RHIGH .GT. TOL) THEN 
INTEGRAL=EVALT(RADIUS) 
IF(ABSCINTEGRAL - 0.5) .LT. TOL) GOTO 75 
IF(INTEGRAL .GT. 0.5) THEN 
RHIGH = RADIUS 
TPHIGH=INTEGRAL 
ELSE 
RLOW=RADIUS 
TPLOW=INTEGRAL 
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ENDIF 
RAD IUS=RLOW+(RHIGH-RLOW)*CRATIO 
IF (ICNT .GE. 150) THEN 
PRINT*, ‘COUNT EXCEEDED ' 
GOTO 75 
ENDIF 
GOTO 60 
ENDIF 
75 CONTINUE 
R=RADIUS 
RETURN 
END 


REAL FUNCTION EVALT(RAD) 
* EVALUATES THE AREA COVERAGE BY TRIVARIATE NORMAL FOR INPUT R 


INCLUDE 'COM DEF' 

INCLUDE 'SEPCOM DEF' 

REAL ANS,LOMWER, UPPER, RAD 

EXTERNAL RAPFC3, RAPLOW, RAPUP 

RADIUS=RAD 

LOWER=0.0 

UPPER=6. 2831853 

CALL DTWODQ(RAPFC3, LOWER, UPPER, RAPLOW,RAPUP, ERRABS, ERRREL, IRULE, 
C ANS, ERREST) 


EVALT=ANS 
RETURN 
END 


REAL FUNCTION RAPLOW(X) 
REAL X 

RAPLOW=0.0 

RETURN 

END 


REAL FUNCTION RAPUP(X) 
REAL X 
RAPUP=3.1415927 
RETURN 

END 


REAL FUNCTION RAPFC3CTHETA,PHI) 
* FUNCTION USED TO FIND SEP 


REAL CP,SP,CT,ST,AA,SA,SAM1,BB,82DA, ERFARG, EXPARG, THETA, PHI 

INCLUDE 'COM DEF' 

INCLUDE 'SEPCOM DEF' 

EXTERNAL ERF 

СР-СО5(РНІ) 

CT=COS(THETA) 

SP=SIN(PHI) 

ST=SIN(THETA) 

AAZ(CXX*SP*SP*CT*CT«CYY*SP*SP*ST*ST*CZZ*CP*CP)/2.0 

SA=SQRT (AA) 

SAM1=1.0/SA 

BB=(CXX*MN(1)*SP*CT+CYY*MN(2)*SP*ST+CZZ*MN(3)*CP)/2.0 

B2DA=BB*BB/AA 

ERFARG=SA*RADIUS+SAM1*BB 

EXPARG=-AA*RADIUS*RADIUS-2.0*BB*RADIUS-B2DA 

IF (EXPARG .GT. -32.0) THEN 
F=ERFCERFARG)*SAM1*(0.5+B2DA) -EXPCEXPARG)* (RADIUS-BB/AA) 
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*SQPIINV 
ELSE 
F=ERF(ERFARG)*SAM1*(0.5+B2DA) 
ENDIF 


ERFARG=SAM1*BB 
EXPARG=-B2DA 
IF (EXPARG .GT. -32.0) THEN 
F=F-ERF(ERFARG)*SAM1*(0.5+B2DA) -EXP(EXPARG)*(BB/AA)*SOPIINV 
ELSE 
F=F-ERF(ERFARG)*SAM1*(0.5+B2DA) 
ENDIF 
EXPARG=B2DA-C 
IF (EXPARG .GT. -32.0) THEN 
F=F*CSANT*EXP(EXPARG)*SP/AA 
IFCF .LT. 0.000000000000001) F=0.0 
ELSE 
F=0.0 
ENDIF 
RAPFC3=F 
RETURN 
END 


SUBROUTINE DTIVAL (DT1,SEPMU,SEPSIG) 


SETS THE VALUES OF THE MATRIX DT1 AND THE VECTORS SEPMU AND 
SEPSIG. 


INCLUDE "СОМ DEF! 
INCLUDE 'WVEC DEF! 
REAL DT1(10,9), SEPMU(3), SEPSIG(6) 
INTEGER I,J 
00 100 1-1,9 
DO 200 7-1,9 
1Е (1 .NE. J) THEN 
DT1(I,J)20.0 
ELSE 
DT1(1,J)=1.0 
ENDIF 
CONT I NUE 
CONTINUE 
DT1(10,1)=W(8)/W(5) 
DT1(10,2)=W(10)/W(5) 
DT1(10,3)=W(6)/(2.0*W(5)) 
DT1(10,4)=(R/(8.0*W(5)))*(4.O*SIGINV(1)*W(9)+2.0*SIGINV(2)*WC11) 


С +SIGINVC(3)*(WC1)-W(2)))-(WCB)/(2.0*W(5)))* 

С (SIGINVC1)*MUC1)*SIGINV(2)*MU(2)* SIGINV(3) *MUC3)) 
DT1C(10,5)2(R/(4.0*W(5)))*C2.0*SIGINV(1)*W(11)*4.0* SIGINV(2)*W( 12) 
C +SIGINVC3)*(WC(3)-WC(4)))-(WC10)/W(5))* 

C (SIGINV(1)*MU(1)*SIGINV(2)*MU(2)*SIGINV(3)*MU(3)) . 
DT1(10,6)=(R/(&.0*W(5)))*(SIGINV(1)*(W(1)-W(2))+SIGINV(2)*(W(3) 
С -W(4))+SIGINV(3)*CW(7)+WC(5)))-(WC(6)/(2.0*W(5)))* 

C (SIGINVC1)*MUC1)*SIGINV(2) *HMUC2) * SIGINV(3)*MU(C3)) 
DT1(10,7)2(R/(8.0*w(5)))*(2.0* SIGINV(2)*W(11)*4 .0* SIGINV(4) 

C *w(12)*SIGINV(5)*(W(3)-W(4)))-CWC105/(2.0*W(5)))* 

С CSIGINV(2)*MUC1)+SIGINV(4)*MU(2)+SIGINV(5)*MU(3)) 


DT1(10,8)=(R/(&.0*W(5)))*(SIGINV(2)*(W(1)-W(2))+SIGINV(&)*(W(3) 


с -W(&))+SIGINV(5)*(W(7)+W(5)))-(W(6)/(2.0*Ww(5)))* 
C 


C 
C 


(SIGINV(2)*MU(1)+SIGINV(&)*MU(2)+SIGINV(5)*MU(3)) 
DT1(10,9)=(R/(8.0*Ww(5)))*(SIGINV(3)*(W(1)-W(2))+SIGINV(5)*(W(3) 
-W(&4))+SIGINV(6)*(W(7)+W(5)))-(W(6)/(&.0*Ww(5)))* 
(SIGINV(3)*MU(1)+SIGINV(5)*MU(2)+SIGINV(6)*MU(3)) 
DO 10 I=1,3 
SEPMU(CI)-DT1(10,1) 


10 CONTINUE 


DO 20 1-1,6 
SEPSIG(1)=DT1(10,1+3) 
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ж 


20 


200 
100 


CONTINUE 
RETURN 
END 


SUBROUTINE DT2VAL(DT2) 


SETS THE VALUES OF THE MATRIX DT2. THE VECTORS A1 AND B1 ARE 
IN THIS SUBROUTINE AND THEN PASSED TO THE SUBROUTINE FIXMAT, 
WHICH ACTUALLY SETS THE VALUES OF DT2. 


INCLUDE 'COM ПЕР" 

INCLUDE 'WVEC DEF' 

REAL DT2(22,10),SIGM(3,3), TEMP(3),ANS(3),B1(3),A1(6) 

INTEGER I,J 

DO 100 1-1,10 
00 200 J=1,10 

IF (I „МЕ. J) THEN 
DT2(1,J)20.0 
ELSE 
DT2(1,J)21.0 
ENDIF 
CONT INUE 

CONT INUE 

SIGM(1,1)=SIGINV(1) 

SIGM(1,2)=SIGINV(2) 

SIGM(1,3)=SIGINV(3) 

SIGM(2, 1)=SIGINV(2) 

SIGM(2,2)=SIGINV(4) 

SIGM(2,3)=SIGINV(5) 

SIGM(3,1)=SIGINV(3) 

SIGM(3,2)=SIGINV(5) 

SIGM(3,3)=SIGINV(6) 

B1(1)=SC(2)/2.0 

B1(2)=SS(3)/4.0 

B1(3)=CC(10) 

А1(12-СС(2)-СС(12) 

A1(2)22.0* (W(3)-CS(17)-CS(2)*CS(11)) 
А1(3)-2.0%(5С(1)-М(9)) 
А1(43:М(1)-СС(17)-СС(2)%СС(12) 
А1(52-55(1)-М(11) 

A1(6)=CC(17) 

CALL FIXMAT(A1,B1,11,DT2,SIGM) 
В1(12-(5С(5)-5С(2))/2.0 
B1(2)=(SS(7)-SS(3))/4.0 
B1(3)=(CC(6)+CC(3))/2.0 
A1C1)=(2.0*CC(5)-CC(9)-CC(2))/4.0 
А1(2)-(2.0%М4(4)-С5(8)-М(3)-2.0%С5(5)%С5(9)%С5(2)2/2.0 
A1(3)=(SC(7)-SC(1))/2.0 
A1(4)2(2.0*W(2) -CC(8)-W(1)-2.0*CC(5)*CC(9)*CC(2))/4.0 
A1(5)=(SS(9)-SS(1))/4.0 
A1(6)=(CC(8)+2.0*W(2)+W(1))/4.0 

CALL FIXMAT(A1,B1, 12,DT2,SIGM) 
B1(1)=SS(3)/4.0 
B1(2)=SS(4)/2.0 
B1(3)=CS(10) 
A1C(1)24(3)-CS(17)-CS(2)*CS(11) 
A1(2)22.0*(W(1)-CC(2)-CC( 17)*CC( 12)) 
А1(3)=55(1) -М(11) 

A1(4)=CS(2)-CS(11) 
А1(5)-2.0%(55(2)-М(12)) 
А1(6)-С5(17) 

CALL FIXMAT(A1,B1, 13,DT2,SIGM) 
B1(1)=(SS(7)-SS(3))/4.0 
B1(2)=(SS(8)-SS(4))/2.0 
B1(3)=(CS(6)+CS(3))/2.0 
A1C1)=(2.0*W(4)-2.0*CS(5)-CS(B)+CS(9)-W(3)+CS(2))/4.0 
A1C2)=(2.0*W(2)-2.0*CC(5)-CC(B)+CCC(9)-WC(1)+CC(2))/2.0 
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А1(3)=(55(9) -55(1))/4.0 
A1(4)=(2.0*CS(5)-CS(9)-CS(2))/4.0 
A1(5)=(SS(10)-SS(2))/2.0 
A1(6)=(CS(8)+2.0*W(4)+W(3))/4.0 

CALL FIXMAT(A1,B1,14,DT2,SIGM) 
B1(1)=W(8) 

B1(2)=W(10) 
B1(3)=W(6)/2.0 
A1C1)=W(9) 

A1(2)=W(11) 
A1(3)=2.0*(W(1)-CC(17)) 
A1(4)=W(12) 
A1(5)=2.0*(W(3)-CS(17)) 
A1(6)=W(5)-SC(9) 

CALL FIXMAT(A1,B1, 15,DT2,SIGM) 
B1(1)=(W(1)-W(2))/2.0 
B1(2)=(W(3)-W(4))/2.0 
B1(3)=(W(7)+W(5))/2.0 
A1(1)=(2.0*SC(2)-SC(5))/4.0 
A1(2)=(2.0*SS(3)-SS(7))/4.0 
A1(3)=(2.0*W(8)+CC(3)-CC(6))/2.0 
A1(4)=(2.0*SS(4)-SS(8))/4.0 
A1(5)=(2.0*W(10)+CS(3)-CS(6))/2.0 
A1(6)=(2.0*W(6)+SC(4))/4.0 

CALL FIXMAT(A1,B1,16,DT2,SIGM) 
B1(1)=(CC(3)-CC(6))/2.0 
B1(2)=(CS(3)-CS(6))/2.0 
B1(3)=(SC(4)+W(6))/2.0 
A1(1)=(2.0*SC(3)-SC(1)-SC(7))/4.0 
A1(2)=(2.0*SS(5)-SS(1)-SS(9))/4.0 
A1(3)=(W(1)-CC(8))/2.0 
A1(4)=(2.0*SS(6)-SS(2)-SS(10))/4.0 
A1(5)=(W(3)-CS(8))/2.0 
A1(6)=(SC(6)+2.0*W(7)+W(5))/4.0 

CALL FIXMAT(A1,B1, 17,DT2, SIGM) 
B1(1)=W(9) 

B1(2)=W(11)/2.0 
B1(3)=w(1)-CC(17) 
A1(1)=SC(12) 
А1(2)-2.0%(55(12)-55(14)) 
A1(3)=(2.0*SC(2)-SC(5))/4.0 
А1(4)-5С(10)-5С(12) 
А1(5)-(2.0%55(3)-55(7))/8.0 
A1(6)=W(8)-SC(10) 

CALL FIXMAT(A1,B1, 18,DT2, SIGM) 

B1(1)=SC(12) 

В1(2)-55(12)-55(14) 

B1(3)=(2.0*SC(2)-SC(5))/8.0 

A1(1)=SC(16) 

А1(2)=(55(16)+2.0*55(15))/4.0 
A1(3)22.0*(CC(2)-2.0*CC(12)*CC(16)) 

А1(0)-55(17)-55(18) 
A1(5)=2.0*(W(3)-CS(2)-2.0*CS(17)+2.0*CS(11)+CS(13)-CS(16)) 
A1(6)=W(9)-SC(15) 

CALL ҒІХМАТ(А1,81,19,072,5106М) 
B1(1)=W(11)/2.0 
B1(2)=W(12) 

B1(3)=W(3)-CS(17) 
A1(1)=SS(12)-SS(14) 
A1(2)=2.0*(SC(10)-SC(12)) 
А1(3)=(2.0*$$(3)-$$(7))/8.0 
A1(4)=SS(14) 
A1(5)=(2.0*SS(4)-SS(8))/4.0 
A1(6)=W(10)-SS(12) 

CALL FIXMAT(A1,B1,20,DT2, SIGM) 
81(1)-(55(13)%55(12))/2.0 
B1(2)-(SC(10)-SC(11))/2.0 
B1(3)=(2.0*SS(3)-SS(7))/8.0 
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ж 


С 


100 


200 


500 


A1(1)=(SS(16)+2.0*SS(15))/4.0 
А1(2)=(2.0*$С(15)-$С(14)-$С(13))/2.0 
A1C3)=CS(1)+W(3)-2.0*CS(18)-2.0*CS(17)+CS(14)+CS(13) 
A1(4)=(2.0*SS(15)-SS(16))/4.0 
А1(59-М(1)-СС(1)-2.0%СС(17)%2.0%СС(119%СС(14)-СС(15) 
A1(6)=W(11)-SS(15) 

CALL FIXMAT(A1,B1,21,D12, SIGM) 
B1(1)=SC(10)-SC(12) 
B1(2)=SS(14) 
B1(3)=(2.0*SS(4)-SS(8))/8.0 
А1(1)-55(17)-55(18) 
A1(2)2(2.0*SS(15)-SS(16))/4.0 
A1(3)22.0* (4(1)-CC(2)-2.0*CC(17)*2.0*CC( 12) *CC( 14) -CC( 162) 
А1(0)-55(18) 
А1(5)-2.0%(С5(2)-2.0%С5(11)%С5(16)) 
А1(6)-М(12)-55(17) 

CALL FIXMAT(A1,B1,22,DT2,SIGM) 

RETURN 

END 


SUBROUTINE FIXMAT(A1,B1,N,0T2,SIGM) 


SETS THE VALUES OF DT2, GIVEN A1, B1, SIGMA INVERSE (SIGM) 
AND THE ROW NUMBER TO BE SET AS INPUTS. 


INCLUDE 'WVEC DEF' 
INCLUDE 'COM DEF' 
REAL DT2(22,10),SIGM(3,3),A1(6),B1(3),A2(6),A3(6),A4(3),MSUM(6) 
,G(6),ANS(3) 

INTEGER N,I 
00 100 1-1,3 

A4CI)2B1CI )*R-MUCID*W(N- 10) 
CONTINUE 
CALL MXTMLT(SIGM,A4,ANS,3,3,1) 
DT2(N, 1)=ANS(1) 
DT2(N,2)=ANS(2) 
DT2(N,3)=ANS(3) 
A2(1)=B1¢(1)*MUC1) 
A2(2)=B1(2)*MUC1)+B1¢1)*MUC2) 
A2(3)zB1(3)*MUC1)*B1C 1)*MUC3) 
А2(4)=В1(2) *М0(2) 
A2(5)=B1(3)*MUC2)+B1(2)*MU(3) 
A2(6)=B1(3)*MU(3) 
DO 200 I=1,6 

MSUM(I)=R*A1(I)-A2(I) 
CONTINUE 
CALL MXTMLT(SIGINV,MSUM,ANS,1,6,1) 
DT2(N,10)=-ANS(1) 
A3(1)=WCN- 10)*(MUC 1 )**2) 
A3(2)22.0*W(N- 10) *MU(1)*MU(2) 
A3(3)22.0*WU(N- 10) *MU( 1) *MU(3) 
A3(4)zW(N- 10) * (MU(2)**2) 
A3(5)22.0*WU(N- 10) *MU(2)*MU(3) 
A3(6)zW(N- 10) * (MU( 35) **2) 
DO 300 1-1,6 

MSUM( I )Z(R**2)*A1(I1)-2.0*R*A2(I)*A3(I) 
CONTINUE 
G(1)=SIGINV(1)*SIGINV(1) 
G(2)=SIGINV(2)*SIGINV(1) 
G(3)=SIGINV(3)*SIGINV(1) 
G(4)=SIGINV(4)*SIGINVC(1)-SIG(6)/DET 
G(5)5SIGINV(5)*SIGINV(1)*SIG(5)/DET 
G(6)=SIGINV(6)*SIGINV(1)-SIG(4)/DET 
CALL MXTMLT(G,MSUM,ANS,1,6,1) 
DT2(N,4)=CANS(1)-W(N-10)*SIGINV(1))/2.0 
G(1)=SIGINV(1)*SIGINVC(2) 
G(2)=SIGINV(2)*SIGINV(2)+S1G(6)/(2.0*DET) 
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G(3)=SIGINV(3)*SIGINV(2)-SIG(5)/(2.0*DET) 
G(4)=SIGINV(4)*SIGINV(2) 
G(5)=SIGINV(5)*SIGINV(2)-SIG(3)/(2.0*DET) 
G(6)=SIGINV(6)*SIGINV(2)+S1G(2)/DET 

CALL MXTMLT(G,MSUM, ANS, 1,6, 1) 
DT2CN,5)=ANS(1)-W(N-10)*SIGINV(2) 
G(1)=SIGINV(1)*SIGINV(3) 
G(2)=SIGINV(2)*SIGINV(3)-SIG(5)/(2.0*DET) 
G(3)5SIGINV(3)*SIGINV(3)*SIG(4)/(2.0*DET) 
G(4)=SIGINV(4)*SIGINV(3)+SIG(3)/DET 
G(S)=SIGINV(S)*SIGINV(3)-SIG(2)/(2.0*DET) 
G(6)=SIGINV(6)*SIGINV(3) 

CALL MXTMLT(G,MSUM,ANS,1,6,1) 
DT2(N,6)=ANS(1)-W(N-10)*SIGINV(3) 
G(1)=SIGINV(1)*SIGINV(4)-SIG(6)/DET 
G(2)=SIGINV(2)*SIGINV(4) 
G(3)=SIGINV(3)*SIGINV(4)+SIG(3)/DET 
G(4)=SIGINV(4)*SIGINV(4) 
G(5)=SIGINV(5)*SIGINV(4) 
G(6)=SIGINV(6)*SIGINV(4)-SIG(1)/DET 

CALL MXTMLT(G,MSUM,ANS, 1,6, 1) 
DT2CN,7)=(ANS(1)-W(N-10)*SIGINV(4))/2.0 
G(1)=SIGINV(1)*SIGINV(5)+SIG(5)/DET 
G(2)=SIGINV(2)*SIGINV(5)-SIG(3)/(2.0*DET) 
G(3)=SIGINV(3)*SIGINV(5) -SIG(2)/(2.0*DET) 
G(4)=SIGINV(4)*SIGINV(5) 
G(S)=SIGINV(S)*SIGINV(S5)+SIG(1)/(2.0*DET) 
G(6)=SIGINV(6)*SIGINV(5) 

CALL MXTMLT(G,MSUM,ANS,1,6,1) 
DT2CN,8)=ANS(1)-W(N-10)*SIGINV(S) 
G(1)=SIGINV(1)*SIGINV(6)-SIG(4)/DET 
G(2)=SIGINV(2)*SIGINV(6)+SIG(2)/DET 
G(3)=SIGINV(3)*SIGINV(6) 
G(4)=SIGINV(4)*SIGINV(6) -SIG(1)/DET 
G(5)=SIGINV(5)*SIGINV(6) 
G(6)=SIGINV(6)*SIGINV(6) 

CALL MXTMLT(G,MSUM, ANS, 1,6, 1) 
DT2(N,9)zC(ANSC1)-W(N- 10) *SIGINV(6))/2.0 
RETURN 

END 


SUBROUTINE DT3VAL(DT3) 
SETS THE VALUES OF THE MATRIX DT3 


INCLUDE “СОМ DEF! 

INCLUDE 'WVEC DEF! 

REAL DT3(9, 22) ,U41,U42,U43,U51,U52,U61,U62,U71,U72,U73, U81, U82, 
с U91,U92,U93 

INTEGER I,J 

DO 100 1=1,3 

DO 200 J=1,22 
DT3(1,J)=0.0 
200 CONTINUE 
100 CONTINUE 

DT3(1,15)2-W(8)/(W(5)**2) 

DT3(1,18)=1.0/W(5) 

DT3(2,15)=-W(10)/(W(5)**2) 

DT3(2,20)=1.0/W(5) 

DT3(3,15)=-0.5*W(6)/(W(5)**2) 

DT3(3,16)=0.5/W(5) 
U41-R*(2.0*SIGINV(1)*W(9)*SIGINV(2)*W(11)*0.5*SIGINV(3)* (W(1) 
C -W(2)))/(4.0*W(5)) 

U42=W(8)/(2.0*W(5)) 
U43-SIGINV(1)*MU(1)*SIGINV(2) *HU( 2) *SIGINV(2) *MU(3) 
U4=U41-U42*U43 
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DT3(4,1)=-SIGINV(1)*U42 

DT3(4,2)=-SIGINV(2)*U42 

DT3(4,3)=-SIGINV(3)*U42 
US1sR*(CSIGINV(1)*W(11)*2.0*SIGINV(2)*W(12)*0.5*SIGINV(3)*(W(3) 
C -W(4)))/(2.0*WC(5)) 

US2=W(10)/W(5) 

U5zU51-U52*U43 

DT3(5,1)2-SIGINV(1)*U52 

DT3(5,2)=-SIGINV(2)*U52 

DT3(5,3)=-SIGINV(3)*U52 
U61=R*(SIGINV(1)*(W(1)-W(2))+SIGINV(2)*(W(3)-W(&))+SIGINV(3) 
C *(W(7)+W(5)))/(&.O0*W(5)) 

U62=W(6)/(2.0*W(5)) 

U6=U61-U62*U43 

DT3(6,1)=-SIGINV(1)*U62 

DT3(6,2)=-SIGINV(2)*U62 

DT3(6,3)=-SIGINV(3)*U62 
U71=R*(SIGINV(2)*W(11)+2.0*SIGINV(4)*W(12)+0.5*SIGINV(5) 

С *(W(3)-W(&)))/(&.0*W(5)) 

U72=W(10)/(2.0*W(5)) 
U73=SIGINV(2)*MUC(1)+SIGINV(4)*MU(2)+SIGINV(5)*MU(3) 
U7=U71-U72*U73 

DT3(7,1)=-SIGINV(2)*U72 

DT3(7,2)=-SIGINV(4)*U72 

DT3(7,3)=-SIGINV(5)*U72 
U81=R*C(SIGINV(2)*C(W(1)-W(2))+SIGINV(4)*(W(3)-W(4))+SIGINV(5) 
С *(W(7)+W(5)))/(&.0O0*W(5)) 

U82=W(6)/(2.0*W(5)) 

U8=U81-U82*U73 

DT3(8,1)=-SIGINV(2)*U82 

DT3(8,2)=-SIGINV(4)*U82 

DT3(8,3)=-SIGINV(5)*U82 
U91=R*(SIGINV(3)*(W(1)-W(2))+SIGINV(S)*(W(3)-W(4))+SIGINVC(6) 
С *(W(7)+W(5)))/(8.0*W(5)) 

U92=W(6)/(4.0*W(5)) 

U93=SIGINV(3)*MUC(1)+SIGINVC5 )*MUC2)+SIGINV(6)*MUC3) 
U9=U91-U92*U93 

DT3(9,1)=-SIGINV(3)*U92 

DT3(9,2)=-SIGINV(5)*U92 

DT3(9,3)=-SIGINV(6)*U92 

TEMP1z((R*W(11)/2.0) -MU(2)*W(8))/(2.0*W(5)) 
TEMP2z(0.25*R*(W(1)-WC2)) -MU(3)*W(8))/(2.0*W(5)) 
TEMP3Z(R*W(9)-MU(1)*W(8))/(2.0*W(5)) 

DT3(4,4)=-SIGINV(1)*U4 
DT3(4,5)=-2.0*SIGINV(2)*U4-(SIG(6)*TEMP1-SIG(5)*TEMP2)/DET 
DT3(4,6)=-2.0*SIGINV(3)*U4-(SIG(4)*TEMP2-SIGC(5S)*TEMP1)/DET 
DT3(4,7)=-SIGINV(4)*U4-(SIG(3)*TEMP2-SIG(6)*TEMP3)/DET 
DT3(4,8)=-2.0*SIGINV(5)*U4-(2.0*SIG(5)*TEMP3-SIG(3)*TEMP1-SIG(2) 
С *TEMP2)/DET 
DT3(4,9)=-SIGINV(6)*U4-(SIG(2)*TEMP1-SIG(4)*TEMP3)/DET 
TEMP1=(R*W( 12) -MUC2)*W( 10) )/WC5) 
TEMP22(0.25*R*(W(3)-W(4))-MUC3)*W(10))/W(5) 
TEMP3=(R*W(11)/2.0-MUC1)*WC10))/W(5) 

DT3(5,4)=-SIGINV(1)*US 
DT3(5,5)=-2.0*SIGINV(2)*U5-(SIG(6)*TEMP1-SIG(S)*TEMP2)/DET 
DT3(5,6)=-2.0*SIGINV(3)*US-(SIG(4)*TEMP2-SIG(5)*TEMP1)/DET 
DT3(5,7)=-SIGINV(4)*US-(SIG(3)*TEMP2-SIG(6)*TEMP3)/DET 
DT3(5,8)2-2.0*SIGINV(5)*U5- (2.0*SIG(5)*TEMP3-SIG(3)*TEMP1-SIG(2) 
C *TEMP2)/DET 
DT3(5,9)=-SIGINV(6)*U5-(SIG(2)*TEMP1-SIG(4)*TEMP3)/DET 
DT3(7,4)=-SIGINV(1)*U7-(SIG(5)*TEMP2-SIG(6)*TEMP1)/(2.0*DET) 
DT3(7,5)=-2.0*SIGINV(2)*U7-(SIG(6)*TEMP3-SIG(3)*TEMP2)/(2.0*DET) 
DT3(7,6)=-2.0*SIGINV(3)*U7-(2.0*SIG(3)*TEMP1-SIG(S)*TEMP3-SIG(2) 
С *TEMP2)/(2.0*DET) 

DT3C7,7)=-SIGINV(4)*U7 
DT3(7,8)=-2.0*SIGINV(5)*U7-(SIG(1)*TEMP2-SIG(3)*TEMP3)/(2.0*DET) 
DT3(7,9)=-SIGINV(6)*U7-(SIG(2)*TEMP3-SIG(1)*TEMP1)/(2.0*DET) 
TEMP 1=(R*(W(3)-W(4))/2.0-MUC2)*WC(6))/(2.0*W(5)) 
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TEMP2=(R*(W(7)+W(5))/2.O-MU(3)*W(6))/(2.0*W(5)) 
TEMP3=(R*(W(1)-W(2))/2.O-MU(1)*W(6))/(2.0*w(5)) 
DT3(6,4)=-SIGINV(1)*U6 
DT3(6,5)=-2.0*SIGINV(2)*U6-(SIGC(6)*TEMP1-SIGC(S)*TEMP2)/DET 
DT3(6,6)=-2.0*SIGINV(3)*U6-(SIG(4)*TEMP2-SIGC(S)*TEMP1)/DET 
DT3(6,7)=-SIGINV(4)*U6-(SIG(3)*TEMP2-SIG(6)*TEMP3)/DET 
DT3(6,8)2-2.0*SIGINV(5)*U6- (2.0*SIG(5) *TEMP3-SIG(3)*TEMP1-SIG(2) 
С *TEMP2)/DET 
DT3(6,9)=-SIGINV(6)*U6-(SIG(2)*TEMP1-SIG(4)*TEMP3)/DET 
DT3(8,4)=-SIGINV(1)*U8-(SIG(S)*TEMP2-SIG(6)*TEMP1)/DET 
DT3(8,5)=-2.0*SIGINV(2)*U8-(SIG(6)*TEMP3-SIG(3)*TEMP2)/DET 
DT3(8,6)2-2.0*SIGINV(3)*UB- (2.0*SIG(3)*TEMP1-SIG(5)*TEMP3-SIG(2) 
C *TEMP2)/DET 
DT3(8,7)=-SIGINV(4)*U8 
DT3(8,8)=-2.0*SIGINV(5)*U8-(SIG(1)*TEMP2-SIG(3)*TEMP3)/DET 
DT3(8,9)=-SIGINV(6)*U8-(SIG(2)*TEMP3-SIG(1)*TEMP1)/DET 
DT3(9,4)=-SIGINV(1)*U9-(SIG(S)*TEMP1-SIG(4)*TEMP2)/(2.0*DET) 
DT3(9,5)2-2.0*SIGINV(2)*U9- (2.0*SIG(2)*TEMP2-SIG(3)*TEMP1-SIG(5) 
C *TEMP3)/(2.0*DET) 
DT3(9,6)=-2.0*SIGINV(3)*U9-(SIG(4)*TEMP3-SIG(2)*TEMP1)/(2.0*DET) 
DT3(9,7)=-SIGINV(4)*U9-(SIG(3)*TEMP3-SIG(1)*TEMP2)/(2.0*DET) 
DT3(9,8)=-2.0*SIGINV(5)*U9-(SIG(1)*TEMP1-SIG(2)*TEMP3)/(2.0*DET) 
DT3(9,9)=-SIGINV(6)*U9 
DT3(4,10)=U41/R 
DT3¢5,10)=U51/R 
DT3(6,10)zU61/R 
DT3(7,10)=U71/R 
DT3(8,10)=081/R 
DT3(9,10)=U91/R 
DO 300 1=4,9 

DO 400 J=11,22 

DT3(1,J)=0.0 
400 CONTINUE 
300 CONTINUE 

DT3(4,11)=(R*SIGINV(3))/(8.0*W(5)) 
DT3(4,12)=-DT3(4,11) 
DT3(4,15)=-U4/W(5) 
DT3(4,18)=-U043/(2.0*W(5)) 
DT3(4,19)=(R*SIGINV(1))/(2.0*W(5S)) 
DT3(4,21)=(R*SIGINV(2))/(4.0*W(5)) 
DT3(5,13)=(R*SIGINV(3))/(4.0*W(5)) 
DT3(5,14)=-DT3(5,13) 
DT3(5,15)2-U5/W(5) 
073(5,20)--043/М(5) 
DT3(5,21)=DT3(4,19) 
DT3(5,22)=(R*SIGINV(2))/W(5) 
DT3(6,11)=(R*SIGINVC1))/(4.0*W(5)) 
DT3(6,12)=-DT3(6,11) 
DT3(6,13)=DT3(4,21) 
DT3(6,14)=-DT3(6,13) 
DT3(6,17)=(R*SIGINV(3))/(4.0*W(5)) 
DT3(6,15)=-U6/W(5)+DT3(6,17) 
DT3(6,16)=-U043/(2.0*W(5)) 
DT3(7,13)=(R*SIGINV(5))/(8.0*W(5)) 
DT3(7,14)=-DT3(7,13) 
DT3(7,15)=-U7/W(5) 
DT3(7,20)=-U073/(2.0*W(5)) 
DT3(7,21)=DT3(4,21) 
DT3(7,22)-(R*SIGINV(4))/(2.0*W(5)) 
DT3(8,11)=DT3(4,21) 
DT3(8,12)=-DT3(4,21) 
DT3(8,13)=DT3(7,22)/2.0 
DT3(8,14)=-DT3(8,13) 
DT3(8,15)=-08/W(5)+DT3(7,13)*2.0 
DT3(8,16)=DT3(7,20) 
DT3(8,17)=DT3(7,13)*2.0 
DT3(9,11)=DT3(4,11) 
DT3(9,12)=-DT3(9,11) 
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20 


40 
30 


ж 


20 
10 


40 


0T3(9, 13)=DT3(7, 13) 

DT3(9, 14)=-DT3(9, 13) 
DT3(9,17)-(R*SIGINV(6))/(8.0*W(5)) 
DT3(9,15)2-U9/W(5)*DT3(9, 17) 
DT3(9,16)2-U93/(4.0*W(5)) 

RETURN 

END 


SUBROUTINE MAKES 


CREATES THE MATRICES S AND ST 


INCLUDE 'PMUCOM DEF' 
INTEGER I,J 
DO 10 1-1,9 
DO 20 J=1,6 
$(1,4)=0.0 
CONTINUE 
CONTINUE 
$(1,1)=1.0 
5(2,2)-0. 
$(3,3) 
$(4,2 
5(5,4 
5(6,5 
S(7,3 
s(8,5 
S(9,6 
00 30 1-1,6 
DO 40 2-1,9 
STCI,J)5S(J, 1) 
CONTINUE 
CONTINUE 
RETURN 
END 


5 
19 
:9 
0 
29 
59 
29 
20 
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SUBROUTINE PMUS(INP,OUTP1,0UTP2) 


COMPUTES THE VALUES OF PMU AND PSIG 


INCLUDE 'PMUCOM DEF! 
INCLUDE 'COM DEF' 
REAL INP(3,3),OUTP1(3,3), PROD(9,9), TEMP(9, 6), OUTP2(6, 6) 
INTEGER I,J 
CALL TENSRCINP, INP, PROD,3,3,3,3) 
CALL MXTMLT(PROD,S,TEMP,9,9,6) 
CALL MXTMLTCST,TEMP,OUTP2,6,9,6) 
DO 10 1=1,6 
DO 20 J=1,6 
OUTP2(I,J)2(2.0*OUTP2(I, J))/(REAL(CNUMTR) *NUMSIM) 
CONT I NUE 
CONTINUE 
DO 30 1=1,3 
DO 40 J=1,3 
OUTP1(1,J)=INPC1,J)/REALCNUMTR) 
CONT I NUE 
CONTINUE 
RETURN 
END 
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SUBROUTINE TENSR(A,B,C,A1,A2,81,82) 
ж COMPUTES THE TENSOR PRODUCT OF A AND B OUPTUT TO C 


INTEGER I,J,K,L,A1,A2,B1,B2 

REAL A(A1,A2),B(B1,B2),C(A1*B1,A2*B2) 

00 10 1-1,А1 

DO 20 J=1,A2 
00 30 К-1,81 
DO 40 L=1,B2 
C(K+(1-1)*B1,L+(J-1)*B2)=A(1,J)*B(K,L) 

40 CONTINUE 
30 CONTINUE 
20 CONTINUE 
10 CONTINUE 

RETURN 

END 


SUBROUTINE BNDEST(SEPM,PMU,PSIG,SEPMU,SEPSIG,CI) 


* CALCULATES THE FIRST AND SECOND ORDER VARIANCE ESTIMATES AND 
* CORRESPONDING CONFIDENCE INTERVAL SIZES 


INTEGER I,J 
REAL TEMP1(3,3), TEMP2(6,6), TEMP3(3,6), TEMP4(6,3),SSQMU(9), 
C SSQS1G(36) , SSQMSG(18) , SSQSGM( 18) , SEPM(9, 9) , PSIGSQ(36, 36), 
C PSIG(6,6),SEPMU(3), SEPSIG(6), TEMP(36) ,PMUSQ(9,9), PMU(3,3), 
C PMUSIG(18, 18) , PSIGMU(18, 18), CI(2) , STPMU(9) , STPSIG(36) 
C , BND,BND1,BND2,BND3 
DO 10 1=1,3 
DO 20 J=1,3 
TEMP1(I,J)sSEPM(I, J) 
20 CONT INUE 
00 30 1-1,6 
TEMP3(1, J)=SEPM(1, J+3) 
30 CONTINUE 
10 CONTINUE 
CALL STRING(TEMP1,SSQMU, 5,3) 
CALL STRING(TEMP3,SSQMSG, 3,6) 
00 40 1-1,6 
DO 50 J=1,3 
TEMP4(1, J)=SEPM(I+3, J) 
50 CONTINUE 
DO 60 J=1,6 
TEMP2(1,J)=SEPM(1+3,J+3) 
60 CONTINUE 
40 CONTINUE 
CALL STRING(TEMP4, SSQSGM, 6, 3) 
CALL STRING(TEMP2, SSQSIG, 6,6) 
CALL MXTMLT(PMU, SEPMU, TEMP, 3,3, 1) 
CALL MXTMLT(SEPMU, TEMP , BND, 1,3, 1) 
CALL MXTMLT(PSIG,SEPSIG,TEMP,6,6, 1) 
CALL MXTMLT(SEPSIG, TEMP, BND1, 1,6, 1) 
BND1=BND1+BND 
CALL TENSR(PMU,PMU, PMUSQ,3,3,3,3) 
CALL MXTMLT(PMUSQ,SSQMU,TEMP,9,9, 1) 
CALL MXTMLT(SSQMU, TEMP, BND, 1,9, 1) 
CALL TENSR(PSIG,PSIG,PSIGSQ,6,6,6,6) 
CALL MXTMLT(PSIGSQ, SSQSIG, TEMP, 36, 36, 1) 
CALL MXTMLT(SSQSIG, TEMP, BND2, 1,36, 1) 
BND2=(BND+BND2)/2.0 
CALL TENSR(PMU,PSIG,PMUSIG,3,3,6,6) 
CALL MXTMLT(PMUSIG, SSQMSG, TEMP, 18, 18, 1) 
CALL MXTMLT(SSQMSG, TEMP, BND, 1, 18, 1) 
CALL TENSR(PSIG, PMU, PSIGMU,6,6, 3,3) 
CALL MXTMLT(PSIGMU, SSQSGM, TEMP, 18,18, 1) 
CALL MXTMLT(SSQSGM, TEMP, BND3, 1,18, 1) 
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20 
10 


BND2=BND2+( (BND+BND3)/2.0) 

CALL STRINGCPMU,STPMU,3,3) 

CALL STRING(PSIG,STPSIG,6,6) 

CALL MXTMLT(SSQMU,STPMU, BND, 1,9, 1) 
CALL MXTMLT(SSQSIG,STPSIG, BND3, 1,36, 1) 
BND2=BND2+(((BND+BND3)**2)/4.0) 
BND=BND1+BND2 
C1(1)=1.96*SQRT(BND1) 
C1(2)=1.96*SQRT (BND) 

RETURN 

END 


SUBROUTINE STRING(A1,A2,N1,N2) 


STRINGS OUT THE MATRIX A1 AS A VECTOR A2 


REAL A1(N1,N2),A2(N1*N2) 

INTEGER I,N1,N2, J 

DO 10 I=1,N1 
DO 20 J=1,N2 

A2((I-1)*N2+J)=A1(I1I,J) 

CONTINUE 

CONT INUE 

RETURN 

END 


REAL FUNCTION HX(X,Y) 
REAL X,Y,SPHER(3), TEMP 
INCLUDE "СОМ DEF! 
SPHER( 1) -R*SIN(X)*COS(Y) -MU(1) 
SPHER(2)sR*SIN(X)*SIN(Y)-MU(2) 
SPHER(3)-R*COS(X)-MU(3) 
TEMPz-0.5*(((SPHER(1)**2)*SIGINV(1))*(2. 0* SPHER( 1)* SPHER(2) 
*SIGINV(2))*(2.0*SPHER(1)*SPHER(3)*SIGINV(3)) 
*((SPHER(2)**2)*SIGINV(4))*(2.0* SPHER(2)*SPHER(3) 
*SIGINV(5) )*((SPHER(3)**2)*SIGINV(6))) 
IF (TEMP .GT. -27.0 ) THEN 
HXsEXP (TEMP ) /DENOM 
ELSE 
НХ=0.0 
ENDIF 
RETURN 
END 


REAL FUNCTION CC1111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC1111=COS(X)*COS(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC1311(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC1311=COS(3.0*X)*COS(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC3111(X,Y) 
EXTERNAL HX 

REAL X,Y 

CC3111=(COS(X)**3)*COS(Y )*HX(X,Y) 
RETURN 

END 
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REAL FUNCTION CS1111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS1111=COS(X)*SINCY )*HX(X, Y) 
RETURN 

END 


REAL FUNCTION CS1311(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS$1311=COS(3.0*X)*SINCY )*HX(X, Y) 
RETURN 

END 


REAL FUNCTION CS3111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS3111=(COS(X)**3)*SIN(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1100(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC1100=SIN(X)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1200(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC1200=SIN(2.0*X)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1300(X,Y) 
EXTERNAL HX 

REAL X,Y 

SC1300=SIN(3.0*X )*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC2111(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC2111=(SIN(X)**2)*COS(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC3121(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC3121=(SIN(X)**3)*(COSC(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS2111(X,Y) 
EXTERNAL HX 

REAL X,Y 
SS2111=(SIN(X)**2)*SIN(Y)*HX(X,Y) 
RETURN 

END 
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REAL FUNCTION SS3112(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS3112=(SIN(X)**3)*SIN(2.0*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS3121(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS3121=(SIN(X)**3)*(SIN(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC1113(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC1113=COS(X)*COS(3.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC1131(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC1131=COS(X)*(COS(Y)**3)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC1211(X, Y) 
EXTERNAL HX 

REAL X,Y 
CC12112C0S(2.0*X)*COS(Y ) * HX(X, Y) 
RETURN 

END 


REAL FUNCTION CC1331(X,Y) 

EXTERNAL HX 

REAL X,Y 
CC13312C0S(3.0*X)*(COS(Y) * * 3) *HX(CX,, Y) 
RETURN 

END 


REAL FUNCTION CC1411(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC14112COS(4.0*X)*COSCY)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC1511(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC1511=Cos(5.0*X)*COS(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC1531(X,Y) 

EXTERNAL HX 

REAL X,Y 
CC1531=COs(5.0*X)*(COSC(Y)**3)I*HXC(X,Y) 
RETURN 

END 
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REAL FUNCTION CC2111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC21112(COS(X)**2)*COS (Y) *HX(X, Y) 
RETURN 

END 


REAL FUNCTION CC3113(X,Y) 

EXTERNAL HX 

REAL X,Y 
CC31132(COS(X)* *2)*COS (2.0* Y) *HXCX, Y) 
RETURN 

END 


REAL FUNCTION CC3131(X,Y) 

EXTERNAL HX 

REAL X,Y 
CC3131-2(COS (X) **3) *(COS(Y) **3) *HX(X, Y) 
RETURN 

END 


REAL FUNCTION CC5111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CC5111=(COS(X)**5)*COSCY)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC5113(X,Y) 

EXTERNAL HX 

REAL X,Y 
CC5113=(COS(X)**5)*COS(3.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CC5131(X,Y) 

EXTERNAL HX 

REAL X,Y 
CC51312(COS(X)**5)*(COS(Y)**3) *HX(X, Y) 
RETURN 

END 


REAL FUNCTION CS1113(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS1113=COS(X)*SIN(3.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS1131(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS1131=COS(X)*(SINCY)**3)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS1211(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS$1211=COS(2.0*X)*SINCY)*HX(X,Y) 
RETURN 

END 
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REAL FUNCTION CS1331(X,Y) 

EXTERNAL HX 

REAL X,Y 
CS1331=COS(3.0*X)*(SIN(Y)**3)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS1411(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS1411=COS(4.O*X)*SIN(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS1511(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS15112COS(5.0*X)*SIN(Y)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION CS1531(X,Y) 

EXTERNAL HX 

REAL X,Y 
CS1531=COS(5 .0*X) *(SINCY)**3)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS2111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS2111=(COS(X)**2)*SIN(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS3113(X,Y) 

EXTERNAL HX 

REAL X,Y 
CS3113=(COS(X)**3)*SIN(3.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS3131(X,Y) 

EXTERNAL HX 

REAL X,Y 
CS3131=(COS(X)**3)*(SINCY)**3)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION CS5111(X,Y) 
EXTERNAL HX 

REAL X,Y 
CS5111=(COS(X)**S5)*SIN(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION CS5113(X,Y) 

EXTERNAL HX 

REAL X,Y 
CS5113=(COS(X)**S)*SIN(3.O*Y)*HX(X,Y) 
RETURN 

END 
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REAL FUNCTION CS5131(X,Y) 

EXTERNAL HX 

REAL X,Y 
CS51312(COS Q0 **5) * (SIN(Y)**3) *HX (X , Y) 
RETURN 

END 


REAL FUNCTION SC1121(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC1121=SIN(X)*(COS(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1221(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC1221=SIN(2.0*X)*(COS(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1321(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC1321=SIN(3.0*X)*(COS(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1400(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC1400=SIN(4.0*X)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1421(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC1421=SIN(4.0*X)*(COS(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1500(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC1500zSIN(5.0*X)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION SC1521(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC1521:SIN(5.0*X)*(COS(Y)**2)* HX(X, Y) 
RETURN 

END 


REAL FUNCTION SC3100(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC3100=(SIN(X)**3)*HX(X,Y) 
RETURN 

END 


85 


REAL FUNCTION SC4111(X,Y) 
EXTERNAL HX 

REAL X,Y 
SC4111=(SIN(X)**4)*COS(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC4113(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC4113=(SIN(X)**4)*COS(3.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC4131(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC4131=(SIN(X)**4)*(COS(Y)**3)*HXC(X,Y) 
RETURN 

END 


REAL FUNCTION SC5112(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC5112-(SIN(X)**5)*COS(2.0*Y)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION SC5114(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC5114=(SIN(X)**5)*COS(4.0*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC5121(X,Y) 

EXTERNAL HX 

REAL X,Y 
SCS121=(SIN(X)**5)*(COS(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC5141(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC5141=(SIN(X)**5)*(COS(Y)**4)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SC1212(X,Y) 

EXTERNAL HX 

REAL X,Y 
$C1212-SIN(2.0*X)*COS(2.0*Y)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION SC1412(X,Y) 

EXTERNAL HX 

REAL X,Y 
SC1412-SIN(4.0*X)*COS(2.0*Y)*HX(X, Y) 
RETURN 

END 
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REAL FUNCTION SS1112(X,Y) 
EXTERNAL HX 

REAL X,Y 
SS1112=SIN(X)*SIN(2.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS1121(X,Y) 
EXTERNAL HX 

REAL X,Y 
SS1121=SIN(X)*(SIN(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS1212(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS1212:SIN(2.0*X)*SIN(2.0*Y)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION SS1221(X,Y) 

EXTERNAL HX 

REAL X,Y 

SS1221-SIN(2.0*X)* (SIN(Y)**2)*HX(X , Y) 
RETURN 

END 


REAL FUNCTION SS1312(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS1312-SIN(2.0*X)*SIN(2.0*Y)*HX(X, Y) 
RETURN 

END 


REAL FUNCTION SS1321(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS1321-SIN(3.0*X)*(SIN(Y)* *2)* HX(X , Y) 
RETURN 

END 


REAL FUNCTION SS1412(X, Y) 

EXTERNAL HX 

REAL X,Y 
$S1412-SIN(4.0*X)*SIN(2.0*Y )*HX(X, Y) 
RETURN 

END 


REAL FUNCTION SS1421(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS1421=SIN(4.0*X)*(SIN(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS1512(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS1512zSIN(5.0*X)*SIN(2.0*Y)*HX(X, Y) 
RETURN 

END 
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REAL FUNCTION SS1521(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS1521=SIN(5.0*X)*(SIN(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS4111(X,Y) 
EXTERNAL HX 

REAL X,Y 
SS4111=(SIN(X)**4)*SIN(Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS4113(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS4113=(SIN(X)**4)*SIN(3.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS4131(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS4131=(SIN(X)**4)*(SIN(Y)**3)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS5112(X,Y) 

EXTERNAL HX 

REAL X,Y 

SS5112=(SIN(X)**5 )*SIN(2.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS5114(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS5114=(SIN(X)**S)*SIN(4.O*Y)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS5121(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS5121=(SIN(X)**5)*(SIN(Y)**2)*HX(X,Y) 
RETURN 

END 


REAL FUNCTION SS5141(X,Y) 

EXTERNAL HX 

REAL X,Y 
SS5141-(SIN(X)**5)*(SINCY) **4) *HX(X, Y) 
RETURN 

END 


COM DEF 
REAL R,SIG(6),SIGINV(6),DET,DENOM,MU(3) ,P1,A,B, ERRABS, ERRREL 
C ,MSIG(G, 3), MUSET(3), SIGSET(3,3), SEPSET 
INTEGER IRULE, ISEED, NUMTR 
COMMON/COM1/R,SIG,SIGINV,DET ,DENOM, MU, ERRABS,, ERRREL , MSIG, 
C MUSET, SIGSET,SEPSET, IRULE, NUMTR, I SEED 
PARAMETER(PI=3. 1415927) 
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WVEC DEF 
REAL W(12),CC(17),CS(18),SC( 19) ,SS( 18) 
COMMON/WCC/W,CC,CS,SC,SS 


PMUCOM DEF 
REAL S(9,6),ST(6,9) ,NUMSIM 
COMMON/COM3/S,ST,NUMSIM 


SEPCOM DEF 

REAL CXX,CYY,CZZ,SQPIINV,QUPITW,CSANT,C,RADIUS,MN(3) 
COMMON/COM2/CXX, CYY,CZ2Z,CSANT,C,RADIUS,MN 

PARAMETER (SQPI INV=0.5641896, QUPI TW=0.05626977) 
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